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PGAC Transactions Policies 


It may be noted that some of the papers in this 
issue were originally conceived over a year ago. 
Four of them were presented in preliminary form 
at WESCON in 1956. This long delay between 
conception and publication was partly due to our 
policy of having potential papers criticized by sev- 
eral reviewers. Their comments were edited and 
returned to the authors with suggestions for revis- 
ing the papers. 

Most of the reviewers have been prompt in com- 
menting on the papers—although a few have never 
responded. However, all of the cooperative review- 
ers of automatic control papers have been par- 
ticularly critical as they are instructed to be. In 
fact, the main purpose of the reviews is to eliminate 
papers which are not well written, or are not con- 
tributions to the control field. The papers may be 
tutorial, they may discuss an application of theory 
to the control system, or they may concern an en- 
tirely new development in the field. In any case, 
the PGAC hopes to establish a reputation for 
publishing papers of high quality. 

Ordinarily, a paper is returned to the author 
only once, but if it is a tutorial paper, or if a major 
change is recommended in a potentially worth- 
while paper, it may be returned to the author more 
than once. This causes a delay in the publication 
of the paper, but if it improves the quality of the 
paper and makes it more useful, we believe that 
the delay is justified. 

In the future, when the list of competent re- 
viewers grows longer and the editorial procedure 
becomes more efficient, the delay between the re- 
ceipt of papers and their publication will be shorter. 

A list of current reviewers is given in this issue. 


_ Anyone who is interested in keeping pace with 


recent progress in his special field by reviewing 
papers will be included in future lists if he will send 
his name to the PGAC Transactions editor and 
state his particular interest. Of course, reviewers, 
as well as any other contributors, are urged to sub- 
mit their own papers for future publication. New 
papers concerning any phase of automatic control 
are always welcome, particularly if they are con- 
cerned with over-all control system concepts. It 
is hoped that many potential papers will be con- 
ceived at local PGAC chapter meetings and re- 
ported to the TRANSACTIONS editor. 

It should be noted that special TRANSACTIONS 
may contain papers which have not been reviewed. 
These TRANSACTIONS will include only those pa- 
pers presented at PGAC sponsored symposia or 
conferences. This is the procedure similar to that 
followed by the IRE editorial department: papers 
in the PROCEEDINGS are reviewed, but those pub- 
lished in the NATIONAL CONVENTION RECORD and 
WESCON ConvVENTION RECORD are not. Parts of 
each of these Records willbe distributed to PGAC 
members each year in addition to the regular 
PGAC TRANSACTIONS, special TRANSACTIONS, 
Newsletters, and Bulletins. 

Tentatively, two regular PGAC TRANSACTIONS 
with reviewed papers are to be issued each year. 
The spring issue will contain a PGAC membership 
directory, a registry of those who identify them- 
selves professionally with the automatic control 
field. Particular months for the publications may 
be scheduled after the flow of incoming papers is 
more consistent. In the future it is hoped that the 
TRANSACTIONS will improve in quality and increase 
in quantity as our Professional Group develops. 

—The Editor. 
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The Issue in Brief 


As in the previous TRANSACTIONS, the papers included here are 
varied in nature and in scope. 

The first paper is an easily read discussion with no mathematics. 
The second concerns an application of the describing function method 
and its relation to the root locus of a missile system. Root loci are 
reviewed and related to pulsed-data systems in the next paper, and 
the theme of pulsed data is continued in the fourth paper where 
multirate switches are used and the mathematics becomes more pre- 
dominate. Mathematics is used extensively in the last two papers 
which are concerned with statistical properties and optimization of 
control systems. 

A brief review of each paper follows. 


Analog or Digital Computer for Process Control?, T. M. Stout— 
Page 3 


This is an excellent introduction to process control and its asso- 
ciated problems which are quite different from those found in the 
frequently discussed machine tool or fire control systems. One of the 
several problems, the choice of an analog or digital computer, is 
treated in detail with respect to use, versatility, flexibility, precision, 
accuracy, speed cost, and reliability. Examples of analog and digital 
computer applications in process control are given. 


Analysis of a Nonlinear Control System for Stabilizing a Missile, 
Leonard Atran—Page 8 


The describing function technique is used to determine graphically 
the relationship among frequency, hysteresis band, and system time 
delays in a missile roll control servoutilizing peripheral, tangentially 
operating jets. 

Relationships governing angular position and rates are found to 
be functions of a steady-state oscillating frequency, control full mag- 
nitude, and missile constants. A comparison is made between the 


root locus and amplitude-phase presentation, and an analog computer_ 


study is discussed. 


Root Locus Method of Pulse Transfer Function for Sampled-Data 
Control Systems, Masahiro Mori—Page 13 


The original Japanese version of this paper was published in 
Automatic Control (Japan) in 1955. Translational difficulties, re- 
views, and rewriting at a distance have delayed the publication of 
the paper in the United States. In the meantime, papers by G. W. 
Johnson, D. P. Lindorff, C. G. A. Nordling, E. I. Jury, and J. G. 
Truxal have discussed root loci of sampled data systems in American 
publications. Professor Mori’s paper has been included to illustrate 


the original contribution to automatic control which was done in 
Japan and to present a useful review, with examples, of the relation- 
ships between the root loci in the s and z planes for continuous and 
sample data systems. 


Input-Output Analysis of Multirate Feedback Systems, George M. 
Kranc—Page 21 


A general analytical technique described in this paper permits 
the extension of z transform methods to sampled data systems con- 
taining synchronized switches which do not operate with the same 
sampling rate. These systems are referred to as multirate sampling 
systems. 

Multirate sampling is used for purely computational purposes, 
for converting data in some flight control systems, or for system com- 
pensation to improve performance. 


Statistical Design and Analysis of Multiply-Instrumented Control 
Systems, Robert M. Stewart—Page 29 


Wiener’s least-square filter theory is applied to common types of 
multiply instrumented control systems. Multiply instrumented con- 
trol systems are those in which more than one sensing instrument is 
used such as a radar or communications system using frequency di- 
versity, inertial guidance systems using accelerometer, velocity, and 
position data, or a multivariable process control system. 


A Time Domain Synthesis for Optimum Extrapolators, Carl W. Steeg, 
Jr.—Page 32 


This paper includes considerable mathematical background 
fundamental to time domain synthesis which can be useful to control 
engineers interested in the design of optimum predictors and the 
study of certain nonlinear systems. Because the techniques are new 
to most control engineers, the mathematics alone may provide a 
challenging and useful study. An example at the end of the paper aids 
in understanding the method. 

The method is a direct solution of the integral equation for the 
optimum predictor in terms of the solution to the integral equation 
for the optimum filter, and it circumvents practical difficulties found 
in the techniques derived by Wiener. A procedure is given for avoid- 
ing the solution of the integral equations in the synthesis of the 
optimum extrapolator. 

The synthesis procedure simplifies extrapolator designs which are 
used when the prediction interval is a nonnegative time function in- 
stead of a prediction interval that is always a fixed future time froma 
continuously varying instant. 


Ce 
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Analog or Digital Computer for Process Control?” 
T. M. STOUTt 


Summary—Although much has been written about the push- 
button or computer-controlled factory, there has been little discussion 
of the relative merits of analog and digital computers for control ap- 
plications. Engineers familiar with computers are aware of the ad- 
vantages and disadvantages of analog and digital computers as tools 
for scientific investigations. They have probably not considered their 
suitability for process control applications. In this paper, we will: 


1) Outline some of the tasks that will probably be assigned to a 
process control computer; 

2) State some of the questions that must be answered before any 
computer is selected; 

3) Review some of the characteristics of analog and digital com- 
puters with special reference to process control requirements; 

4) Describe some systems already installed or about to be in- 
stalled, in which computers are used. 


Following this discussion, it will be evident, as might be antici- 
pated in advance, that each type of computer has characteristics 
which make it preferable for particular kinds of applications. Analog 
computers are fast, simple, and inexpensive for relatively small con- 

_trol jobs, while digital computers have accuracy, versatility, and 
flexibility which adapt them to more complex control jobs. 


INTRODUCTION 
Fields of Application 


HE kind of application and the role assigned to 
Thee computer make quite a difference in the char- 

acteristics that the control computer should have. 
The applications that we have in mind, as we approach 
the comparison of analog and digital computers, are 
found in the chemical, petroleum, and related industries 
such as metals, drugs, cement, soap, and paper manu- 
facturing. 

The common characteristic of these industries is 
their use of chemical and physical processes on essen- 
tially continuous streams of raw materials. 

This paper does not discuss the application of com- 

puters or computer-like equipment for machine tool 

control or operation of automatic assembly machines. 
These applications seem to involve a different set of 
design considerations and are receiving attention else- 
where. 


The Role of the Computer 


For plants of the type we are considering, we can 
visualize a system having the general form shown in 
Fig. 1. This picture, used in various forms by many 
authors, shows a process of some kind which takes in 
raw materials, chemicals, water and energy in various 
forms (such as electricity, steam, gas, or fuel oil). By 
‘a combination of steps, a number of products or by- 
‘products are made, and some waste is created. 


-* Manuscript received by the PGAC, June 3, 1957. 
- ¢ The Ramo-Wooldridge Corp., Los Angeles, Calif. 


gy 4 


PRODUCTS 


RAW MATERIALS 


PROCESSING 
UNIT 


ed 
TRANSDUCERS [ 
oe corer bee 
COMPUTER 


Fig. 1—Schematic diagram of a computer control system. 
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Control of a plant like this is now exercised by meas- 
uring flow rates, pressures, temperatures, levels, and 
other quantities which characterize the state of the 
process. This information, in combination with composi- 
tion data from laboratory analyses or on-line analyzers, 
is used to govern the operation of individual controllers 
at various points in the process. These controllers main- 
tain the flow rates, pressures, temperatures, levels, and 
so on at values which are expected to produce the de- 
sired product in the best possible way. These values, 
or “set points,” are established initially by the process 
designers and are modified as time goes on by the engi- 
neers and operators in charge of the process. 

In the computer-controlled process of the future, the 
various measurements will be supplied to a computer 
through appropriate amplifiers and converters. Calcu- 
lations and decisions leading to selection of the optimum 
set points will be made by the computer, and the results 
of its efforts will be transmitted back to the process. 
The operator will receive information about the process 
from the computer and can modify the operation of the 
computer or, in an emergency, he can exercise control 
as he does today by looking at recorders and manually 
changing the set points of the controllers. 


Jobs for the Computer 


To make things more specific, we may look at the 
hypothetical process shown in Fig. 2. In this process, 
a feed material is treated with a reagent (or two feed 
materials are combined) at elevated pressure and tem- 
perature in a vessel where a chemical reaction takes 
place. If the reaction is exothermic, a coolant will be 
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Fig. 2—A hypothetical process showing interconnection of 
sensing devices, computer, and controllers. 


required to control the temperature; if not, heat may 
have to be supplied to keep the reaction going. 

In a system of this kind, the computer could havea 
variety of jobs: 

1) Data Interpretation: The signals from flow trans- 
mitters, thermocouples, and in fact, from any instru- 
ment, are not direct indications of the quantity being 
measured. Scale factors and calibration curves must be 
applied to obtain the desired information. At a mini- 
mum, the same treatment must be applied to the output 
of composition analyzers; in some cases, simultaneous 
equations must be solved to uncover the desired infor- 
mation. 

2) Mixing or Blending Calculations: The flow rate 
of the reagent may require continuous adjustment in 
order to maintain a fixed ratio of reagent to the active 
component of the feed material in the face of variations 
in the feed rate or composition. Or, where the effective- 
ness of the reagent does not increase in direct proportion 
to the amount used, the optimum ratio might be a func- 
tion of fluctuating material and product costs. 

3) Control of Reaction Conditions: Because the prod- 
uct yield may drop off at very high feed rates, it may 
sometimes be advantageous to throw away part of the 
available feed. Raising the pressure by closing the valve 
at the outlet of the reactor might raise the yield while 
cutting down the feed rate. Raising the temperature 
might speed up the reaction but cut down the equilibri- 
um yield. Where such optimum conditions result from 
the combined effects of technical and economic factors, 
the computer might keep the process operating con- 
tinuously at the optimum points despite fluctuations in 
material feed rates and compositions, costs, environ- 
mental and equipment conditions, and so on. 

In the case of batch operations, the computer might 
be used to enforce a rigid schedule of operations. More 
to the point, however, the computer might control the 
process in response to the actual progress of the reac- 
tion. In similar fashion, the computer might assist in 
starting and shutting down a continuous process. 
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4) Data Processing: The computer in a control sys- 
tem of this kind can be expected to perform some data 
processing or logging functions. Flow can be totalized 
or averaged for the accounting department, for exam- 
ple, and various computed quantities as well as basic 
process data can be made available to the process 
engineers. 

5) Diagnostic and Alarm Functions: The computer 
should be designed to check the operation of the process, 
the measuring instruments, and the control system it- 
self, with suitable alarm signals provided to the operator 
in the case of equipment malfunctions or breakdowns. 


SYSTEM DESIGN 
System Analysts 


The example just discussed gives some idea of the 
kinds of things that can be done with a computer in 
process control. Exactly what can and, equally impor- 
tant, should be done and how to do it, must be estab- 
lished by a thorough study of the problem at hand. 
This study must answer such questions as: 

What are the objectives of the process? 

What are the raw materials and products involved? 
What are their characteristics, flow rates, and 
dollar values? What are the variations in these 
properties? 

What quantities can be measured? Controlled? 

What are the relations between the important quanti- 
ties? 

Where are the equipment bottlenecks? 

What are the dynamic characteristics of the process 
and its equipment? 

Can the quantities necessary for improved control 
be measured with the requisite speed and accuracy? 

What kind of computer should be used? How should 
it be used? 

After answering all of these questions, then it is neces- 
sary to draw up detailed specifications of the system 
components, specify their interconnections, outline 
operating procedures for regular and emergency con- 
ditions, and attend to a myriad of other details leading 
to the installation, test, and operation of the complete 
system. 


Choice of Computer 


The choice of an analog or digital computer, the 
subject of this paper, is therefore only one of the prob- 
lems facing the control system designer. It would be a 
mistake to consider that there is a single answer to the 
question, “Analog or Digital Computers for Process 
Control?” and that this discussion will settle the prob- 
lem for all time. 

There is a children’s book called “Are Dogs Better 
Than Cats?” which might provide the theme for the 
remainder of this paper. In the book, a boy says that. 
cats can’t swim, and a girl says dogs can’t climb trees. 
After several pages of this sort of argument, they finally | 
agree that “Dogs are better for dog things, and cats are. 
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better for cat things.” Much the same conclusion should 
be reached on the question of analog vs digital com- 
puters! 


Use of Computer 


The manner in which the computer will be used has 
an important bearing on its specifications. If the com- 
puter supervises conventional minor-loop controllers, 
as shown in Fig. 3(a), its tasks and requirements are 
quite different from those encountered if the computer 
acts directly upon the valves which govern the flow of 
materials and utilities, as shown in Fig. 3(b). In the first 
case, the instantaneous regulation of the process con- 
ditions is left to conventional feedback controllers of 
the kind now used for this purpose; the computer per- 
forms a supervisory function similar to that performed 
by the operator and process engineers. In the event of 
computer failure, control can return to the present man- 
ual mode. In the second case, the computer must be 
fast enough to keep pace with the disturbances affecting 
the process and reliable enough to remain on the job 
between scheduled maintenance periods. 

Although the supervisory type of computer control 
leaves something to be desired in the way of equipment 
economy, it represents the probable course of develop- 
ment. We see little likelihood that the process industries 
will authorize the removal of conventional controllers 
and the substitution of presently untried control com- 
puters. This reading of the crystal ball affects our rank- 
ings of analog and digital computers in the discussion 
that follows. People with a different view of the future 
will perhaps arrive at a different ranking. 


COMPUTER CHARACTERISTICS 


In choosing between analog and digital computers, 
a wide variety of characteristics should be considered.} 
The following list contains the most important of these 
characteristics, with some comments on the relative 
standing of analog and digital computers. Some of the 
following remarks may be controversial, but others are, 
perhaps, pretty generally accepted. 


Versatility 


Digital computers can, in principle, perform any 
kind of calculation that can be specified by a sequence 
of arithmetic operations and yes-or-no decisions. Al- 
though we may concede that ingenious designers could 
develop an analog computer to perform any job that 
might be named, analog computers are best adapted 
to calculations involving the basic arithmetic and calcu- 
lus operations of addition, subtraction, multiplication 
by constants, integration, and single-variable function 
generation. Multiplication of variables, generation of 
functions of several variables, and other complex opera- 
tions can be performed, but less satisfactorily. 


1M. Rubinoff, “Analog vs digital computers,” Proc. IRE, vol. 
41, pp. 1254-1262; October, 1953. 
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Fig. 3—Methods for connecting the control computer to the process: 
(a) supervising conventional controllers, (b) as substitute for 
conventional controllers. 


Flexibility 


Modification of the control functions performed by 
the digital computer requires only a change in the pro- 
gram. With an analog computer, such modifications 
may involve only changes of the interconnections be- 
tween the elements, or they may mean large-scale sub- 
stitution or addition of new computing elements. 


Precision 


The term “precision” denotes “sharpness of defini- 
tion” and is related to the number of figures used to 
represent physical quantities in numerical form. Al- 
though the precision of digital computations can be 
extended at will, analog computer results can generally 
be read to only three or four figures. Either type is 
probably adequate for control use on this score. 


Accuracy 


The term “accuracy” denotes “correctness” or “con- 
formity with truth.” In a digital computer, inaccuracies 
arise from the necessity for approximate calculation 
schemes; in an analog computer, from inaccuracies or 
drifts in the calibration of its components. The digital 
computer can be made to give results as accurate as the 
original data, but additional errors will inevitably be 
introduced by analog computation. 


Speed 


The speeds of both types of computers are probably 
adequate for the kind of supervisory control of present 
processes described above. The digital computer is 
probably better equipped for very slow operations such 
as long-term integration or data storage, while the ana- 
log computer may have an edge for control of very high- 
speed processes. 

Part of the speed advantage of analog computers 
stems from simultaneous or parallel operation of the 
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basic computing elements. If sequential operation is 
used in order to reduce the amount of equipment, either 
to cut costs or to improve reliability, this advantage 
tends to disappear. 


Input-Output 


Since most of the sensing devices used in present proc- 
esses have analog outputs, the analog computer can be 
connected to the process without use of analog-digital 
converters. This advantage is not so significant if data 
will be converted anyway for digital logging and will 
diminish as digital sensing devices become more com- 
mon. 


Amount of Equipment 


With analog computers, the amount of equipment 
needed is roughly proportional to the amount of com- 
puting to be done. With digital computers, on the other 
hand, a significant amount of equipment is needed to 
perform the first calculation, but the amount of equip- 
ment does not increase as rapidly with the amount of 
computing. 


Cost 


Because the cost depends heavily on the amount of 
equipment needed, it can be expected that analog 
computers will be less expensive when only a little 
computing must be done, and that the digital computer 


will have an edge for jobs requiring extensive computing... 


However, it is possible that analog computers will re- 
quire more tailoring to fit particular jobs, so that engi- 
neering costs associated with analog computer installa- 
tions might be higher. 

In discussing costs, we should be careful to compare 
like entities. The usual analog computer is not equipped 
for alarm and logging functions, and any fair cost com- 
parison should include the cost of the equipment for per- 
forming these functions as well as the strictly computa- 
tional part of the control job. 


Reliability 


Designers of either analog or digital computers face a 
difficult task in providing the kind of reliability desired 
by the process industries. One process engineer of my 
acquaintance has said that 10,000 hours ofuninterrupted 
service should be the goal; this time amounts to almost 
14 months on a 24-hour-a-day basis! Actually, most 
processing units are shut down at regular or irregular 
intervals, and scheduled maintenance could be per- 
formed much more frequently than once every 14 
months. 

Reliability is an inverse function of the amount of 
equipment involved in the system. Small analog com- 
puters may therefore offer the most reliable solution for 
simple control jobs, while digital computers may be 
preferable to an equivalent assemblage of analog com- 
ponents for large and complicated control jobs. 

By attention to details, it should be possible to design 
either type of computer with the necessary reliability. 
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Both types will be subject to unpredictable complete 
failures; the analog computer may be more susceptible 
to unsuspected gradual deterioration which leads to 
faulty operation but not to complete failure. 


. 


APPLICATIONS 


Both analog and digital computers have been or will 
be used in the following control applications. We do 
not have access to the reasoning that went into selec- 
tion of the computer for these applications, but we can 
describe some features of the systems. 


Electric Power Systems 


Several utilities are using analog computers to figure 
continuously the best way to supply electricity to their 
customers.2-4 The problem here is to insure that the 
total load is divided between the available generating 
stations in such a way that the total cost to the utility 
is a minimum. This problem is solved by a continuous 
transfer of load from high-cost stations to low-cost 
stations until the incremental costs of delivered power 
are equal. The equipment required includes means for 
measuring the power output of each station, generating 
signals which represent the corresponding incremental 
costs, and making adjustments to bring these costs 
into agreement. 

In this application, the optimizing criterion is relative- 
ly simple and well understood. Alarm and recording 
functions are already handled in other ways. Analog 
computers appear to be an acceptable choice. 


Batch Hydrogenation 


At the Case Institute of Technology in Cleveland, re- 
search has been conducted on computer control of a 
batch hydrogenation process.* Here infrared and refrac- 
tive index instruments and an analog computer were used 
to direct the process in such a manner that the desired 
product was made in a minimum time. 

Solution of differential equations, a natural assign- 
ment for an analog computer, played a prominent part 
in this control program. Nevertheless, in concluding their 
report, the research workers on this project expressed a 
desire for several improvements in their equipment in- 
cluding: 

1) An algebraic equation solver to process informa- 

tion from the two analysis instruments. 

2) Magnetic drum storage for smoothing composition 

data. 

3) A magnetic drum storage and computing system 

for continuous determination of process param- 
eters. 


?R. H. Travers, “Automatic economic dispatching and load 
control—Ohio Edison system,” Trans. AIEE, vol. 76; January, 1957. 


’K. N. Burnett, D. W. Halfhill, and B. R. Shepard, “A new 
automatic dispatching system for electric power systems,” Trans. 


AIEE, vol. 75; July, 1956. 


4E. J. Kompass, “The ‘early bird’ goes automatic,” Control Eng., : 


vol. 3, pp. 77-83; December, 1956. 
® “Process Automation—Report 1, 1954-1956,” Case Institute 
of Technology, Cleveland, Ohio; September, 1956. 
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4) A time-program controller for automatic start-up 
and shut-down operations. 
These functions, as well as those already included, 
could be performed by a general-purpose digital com- 
puter designed for process control applications. 


Chemical Plants 


Early in 1956, Burroughs and du Pont cooperated in 
an experiment of considerable interest.6 Measurements 
taken in a chemical plant at Niagara Falls were con- 
verted to pulse form, transmitted by telephone wire to 
Philadelphia, and processed in a digital computer. 
Yields and material balance data were sent back to 
Niagara Falls for use by the plant operators. 

Although the computer was not connected for direct 
control of the process, a number of useful things were 
learned from this experiment. Two advantages of digital 
systems were highlighted: the ability of digital data to 
survive long-distance transmission, and the ability of a 
digital computer to check the validity of data and com- 
puted results. 

At a meeting on the subject of computer applications 
in industrial control, Faulkner of Imperial Chemi- 
cal Industries, Ltd., of England spoke about the prob- 
lems that would be involved in controlling an ammonia 
plant.’ He quoted a conversation in which the plant man- 
ager asks the operator to drop the gas composition to 

3 per cent, starting about 4:30 A.M., in order to prepare 
for a maintenance job the next day. He then says: 

“The shift manager is exercising judgment, that is to 
say he is guessing on the basis of past experience. The 
fact is, however, that all the data to enable him to make 
a precise observation are known. It is too much to ex- 
pect him to do that sort of job at 2 o’clock in the morn- 
ing. I know that a small analogue black box could be 
designed for this purpose, but numbers of black-box 
slide-rules are not, in my opinion, the way to do the 
job, because this is only one example of many things of 
that sort which have to be done. A digital computer 

- would do all these various jobs on demand and many 
others like them, and no reasonable analogue computer 
could cope with the mass of information which goes to 
the making of the material efficiency of an ammonia 
plant.” 

Careful reading of his remarks indicates that he is 
thinking here about a computer to assist the manager 
and not to exercise direct control. Later, he talks about 
a chemical process and advocates use of a digital com- 
puter for getting maximum yield from a parallel bank 
of reactors, each of which contains catalyst in a different 
state of activity. 


Petroleum Refining 
Esso Standard has announced that an automatic 


6 E, W. James, J. Johnson, E. W. Yetter, and M. A. Martin, 
“How a systems team engineered a plant computer test,” Control 
Eng., vol. 3, pp. 83-86; November, 1956. _ ice 

7J. J. Faulkner, “Introduction to session on the application of 
peste! computers in industrial control,” Proc. IEE, vol. 103, part 
, pp. 98-99; 1956. ee 
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data logger-computer will be installed on a fluid cata- 
lytic cracker in a new refinery in Cuba. One-hundred- 
eleven process variables will be measured, converted to 
digital form, and logged. Some of these quantities, 
stored in analog form, will be supplied to a programmed 
general-purpose analog computer which will sequen- 
tially determine the following “operating guides”: 


1) Carbon burning rate, 

2) Catalyst circulation rate, 

3) Ratio of catalyst flow to oil flow, 

4) Ratio of feed to reactor catalyst holdup, 
5) Conversion, 

6) Per cent weight of hydrogen in coke, 

7) Ratio of coke to total feed, 

8) Heat duty of fractionator reflux (top), 
9) Heat duty of fractionator reflux (middle), 
10) Regenerator superficial velocity, 
11) Reactor superficial velocity. 


This information will be supplied to the operator for 
his guidance in running the unit; it is not planned to 
use these computed quantit es for direct control. The 
computations will be made hourly but can be repeated 
more frequently on demand. 

Since this system is not yet in operation, it is too 
early to tell whether the proposed combination of opera- 
tional amplifiers, chopper-stabilized electronic multi- 
plier-divider modules, relays, and stepping switches will 
perform as well as expected. If it should be decided 
later to use this equipment for direct control, it may be 
hoped that the selected quantities turn out to be useful 
in control, that the available computing speeds will be 
adequate for control, and that program modifications 
can be made without extensive equipment changes. 


CONCLUSION 


It is noteworthy perhaps that at least one manu- 
facturer of analog computers for engineering use now 
is providing digital input-output by means of punched 
tape or an electric typewriter. This development, con- 
sidered in combination with the fact that any digital 
control computer must work between analog sensing 
devices and analog controllers, suggests that the com- 
parison should be one of degree, not kind. We should 
be asking perhaps, “How much analog and how much 
digital equipment?” not “Analog or digital equip- 
ment?” since all control systems will be a combination 
of both. 

Both analog and digital computers will be used for 
process control, the choice depending on the require- 
ments of the particular job. The digital computer ap- 
pears well-adapted to the changing computational and 
logging functions that will be encountered, for example, 
in large and complicated petroleum or chemical proc- 
esses. Analog computers have found and will continue 
to find application in a variety of simpler control prob- 


lems. 


8 “Will the logger-computer optimize refining?” Control Eng., vol. 
4, pp. 23-29; March, 1957. 
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Analysis of a Nonlinear Control System for 
Stabilizing a Missile’ 


LEONARD ATRANT 


Summary—An autopilot with attitude and rate feedback, repre- 
sentative system lags, and a two-way relay servo with inherent 
hysteresis is considered for roll control of a missile with peripheral, 
tangentially operating jets. 

This type of control system is shown to produce a steady-state 
oscillation. Missile dynamics in the presence of this hunting are de- 
veloped and the relationships governing angular position and rates 
are found to be functions of the oscillation frequency, control force 
magnitude, and missile constants (geometry and weight). 

The describing function technique is utilized to determine graph- 
ically the relationship among frequency, hysteresis band, and system 
time delays. A comparison is made between the root locus and am- 
plitude-phase presentation. An analog computer study of system be- 
havior is presented to illustrate the agreement between the analysis 
and system performance. 


INTRODUCTION 


NE of the foremost considerations in the design 

() of a missile control system is that it be as simple 

and inexpensive as possible. This can be achieved 

by using only those components which are absolutely 

necessary to perform the control task adequately. 

Gyros, resolvers, amplifiers, filter networks relays, and 
actuator hardware should be kept to a minimum. 

This paper shall restrict itself to the discussion of 
roll stabilization by means of a relatively simple and 
reliable autopilot system designed to provide attitude 
control, whereas the missile maneuver commands in 
this case are applied to the pitch and yaw channels. 

Although in aircraft operation almost any continuous 
oscillation is intolerable, in an unmanned missile it is 
possible that a steady-state oscillation confined to 
reasonable limits is entirely satisfactory. The control 
system under discussion consists in part of a relay servo 
driving a jet nozzle valve in such a manner as to pro- 
duce a continuous dither. The nozzles themselves are 
located on the periphery of a plane normal to the main 
thrust line. Control of the jets results in a fixed torque 
in either of opposite directions. 

The autopilot signal required to roll stabilize is derived 
from an attitude gyro; however, a derivative signal is 
necessary to augment the negligible aerodynamic damp- 
ing of the missile. 


MiIssILE Rott DyNAmics 


Let us assume a general missile configuration as shown 
in Fig. 1. 

The roll moment equation, taken about the body axis 
is: 


* Manuscript received by the PGAC, September 9, 1956. Pre- 
sented at WESCON, Los Angeles, Calif., August, 1956. 
+ Air Arm Div., Westinghouse Electric Corp., Baltimore, Md. 


Fig. 1—Missile coordinate system. 


T,P + Uz — Iy)QR — InR + PQ) 


= >> L (summation of roll moments). (1) 
Because of symmetry, we may assume 
fy ep ant, =. 0, 
Eq. (1) in simplified form becomes: 
P = lif; + Inf (2) 


where 
1;=L;/Iz (control moment effectiveness), 
lm =Lm/Iz (misalignment torque factor), 
f; =force produced by control jets, 
fm=net unbalanced force due to possible thrust mis- 
alignment, 
L;=control force moment arm, 
L»=effective moment arm upon which unbalanced 
forces act. 
With no unbalanced torques acting on the system, 
fm =0; therefore, upon integrating (1), we arrive at the 
following expression for roll rate: 


P = 1,f;|t — t switching] = ¢ (missile roll rate). 


(3) 


During any “on-time” 


f; = + F (pounds) = constant. 


The on-time for each unidirectional blast is 7/w 
where w is the frequency of oscillation in rad/sec. The 
peak roll rate is, therefore: 


— 901; 


PF 
= ——— (deg/sec). (4) 


max 
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Fig. 2—(a) Waveforms of force and acceleration. (b) Waveforms 
of roll position and rate. 
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Fig. 3—Plot of maximum roll angle and roll rate vs frequency of 
acceleration. 


Further integration leads to the following expression 
for the peak roll angle: 


—711,;F 


ea ok (deg). (5) 


dmax = 
The waveforms of the system under steady-state 
conditions are illustrated in Fig. 2(a) and 2(b). 
Plots of @max and ¢max vs mode frequency are shown 
in Fig. 3 for values of /;=—0.58 and F=100 pounds. 
Superposition of an unbalanced torque requires the 
development of a counter moment of the same average 
value, thereby necessitating that the symmetry of the 
control force vs time curve be destroyed. This is illus- 
trated in the REAC recordings of Fig. 4. 
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Ly SYSTEM. LAG = 0.08 nace 
HET CONTROL FORCE <l00 Ibs ~~ 
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Fig. 4—Relay servo controlling missile in roll (externally 
applied torque of 49 pounds-feet). 


SERVOSYSTEM 


The input-output waveforms of the relay-actuator 
combination are illustrated in Fig. 5, next page. 

The autopilot control system is shown in block dia- 
gram form in Fig. 6. 

In an effort to determine the behavior of the system 
in the presence of a highly nonlinear servo, the system 
transfer functions are combined and the characteristic 
equation rearranged to the following expression for neu- 
tral stability: 


—1 13 1;B(ris +- 1) 
N(x) Sirgen Gaeen ty 


where 7;=A/B. 

The describing function N(x) has been derived in 
the Appendix; in addition, both sides of (6) are plotted 
in Fig. 7. The crossover point of the describing func- 
tion curve and system transfer function determines 
both the frequency and amplitude of the resultant 
steady-state oscillation. 


(6) 


AUTOPILOT PARAMETERS 


It has been shown that for the nonlinear system 
under discussion the magnitude of the roll angle excur- 
sions is a function of the control force and the frequency 
of oscillation. In turn, the frequency is a function of the 
hysteresis band, system time lags, and autopilot gain. 
The interrelation will now be explored. 

A steady-state torque unbalance, due in part to a 
misalignment of the thrust axis, will compel the missile. 
to assume a steady-state angular error. However, the 
system will continue to oscillate in accordance with the 


1E. C. Johnson, “Sinusoidal analysis of feedback control systems 
containing nonlinear elements,” Tvans. AIEE, vol. 71, part 2, pp. 
169-181; July, 1952. 
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Fig. 5—Control system waveforms. 
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Fig. 6—Control system block diagram. 


requirements imposed by (6), and the bias from its origi- 
nal zero steady-state roll angle may be determined in the 
following manner. 

Assume for the moment that the roll oscillation caused 
by system nonlinearity is sufficiently removed in fre- 
quency from the actual roll control mode corner fre- 
quency for complete decoupling to exist. That is, the 
response to a roll command input is not influenced by 
the higher frequency dither. Under these conditions we 
may write that 


and average 
Bisse — BOaver 
Substituting the above into (2) we have: 


IL, 
L; 


RVG gaa ane = mi rad). il 
o " Jm (rad) (7) 

Thus if the maximum allowable roll angle is specified 
together with the magnitude of the unbalanced torques 
likely to be encountered, the required autopilot gain 


becomes: 


Linh 


JPAV 


B= 


(pounds/rad). 


To counterbalance the moment Lmfm requires a force 
capability of 
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AMPLITUDE AND PHASE OF 


=I 
N(X 


vs. w 
) 


L INTERSECTION OF SERVO 
DESCRIBING FUNCTION (A) 
CURVE WITH SYSTEM 
TRANSFER FUNCTION (B) 
EC DETERMINES FREQUENCY 
AND MAGNITUDE OF 
OSCILLATION 


/SEC 


(db) 


Amplitude 


ee | 2 


Fig. 7—Amplitude and phase of servo describing function and 
system transfer function. 


-30 aul: 


De rn 
Fo= jae (pounds). 


“) 


The roll mode, regarded as linear and of the form 
s?—Al;s — Bl; =0, has a damping factor 


sae Ye 


2Wn, 


where 


Wn = W—Bi;. 


Thus for a specified ¢, the value of A that is neces- 
sary is 


Deas 
‘arms f ‘ 

—l; 
DISCUSSION 


A missile control system was selected for analysis 
and the following values were assigned to the various 
parameters. 


1957, 


jes 
——— ie} 
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Fig. 8—Root locus plot of representative control system showing 
effect of nonlinearity. A =linear system, B=nonlinear system. 
72=0.05 sec, F=10h. 


1 
pound — sec? 
F = 100 pounds, 
A = 42 pounds/rad/sec, 
B = 420 pounds/rad, 
Tz = 0, 0.05, 0.1 second, 
T3 = 0.01 second for rz = 0.05, 
h = 0.1F,-0.5F. 


a Os 


L= 


Time lags rz and 73 were considered as representative 
delays for this type system. 

Autopilot and missile characteristics, as illustrated 
by Figs. 7-9, result in the following observations. 

To limit the amplitude of the roll angle oscillation to 
+3.5° requires the missile to oscillate at approximately 
5.5 cps. Thus an autopilot lag of 0.05 second and a 
hysteresis/force ratio of 0.1 will satisfy this condition. 
An increase in either time delay or hysteresis results 
in a lowering of the natural frequency and a consequent 
increase in the signal amplitude level, a factor of im- 
portance for a system operating close to saturation. 

The root locus plot,? while providing the requisite 


2 W. R. Evans, “Graphical analysis of control systems,” Trans. 


AIEE, vol. 67, p. 85; January, 1948. 
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Fig. 9—Relay servo controlling missile in roll. 


stability information, proves to be somewhat more diffi- 
cult to construct as a result of this particular non- 
linearity. The actual procedure consisted of obtaining 
a family of loci for a series of describing function phase 
angles, then plotting on each “phase” locus the gain 
corresponding to the describing function amplitude. 
To analyze the effects of additional lags demands a 
complete new family of loci, whereas the amplitude- 
phase presentation requires but a modified plot of the 
system transfer function. 

A computer (REAC) study resulted in reasonable 
agreement between the analysis and actual system 
performance, the frequency and amplitude checking 
to within 10 per cent thereby justifying the assumption 
regarding the filtering of higher harmonics. 

For example, in Figs. 7 and 8, a resonant point 
exists at w=36 rad/sec when 72=0.05 sec and F=10h 
for the system under discussion. Fig. 9 illustrates the 
results of a REAC study with agreement as to frequency 
of oscillation; the theoretical amplitudes may be found 
in Fig. 3. 


CONCLUSION 


1) Unbalanced torques should be minimized to re- 
duce the control force. 


2) The frequency of the roll oscillation must be main- 
tained at a reasonably high value so as to restrict 
the magnitude of the roll oscillation and to reduce 
coupling with the roll control mode. 


3) Hysteresis and time lags due to nonideal com- 
ponents must be minimized within the autopilot 
in order for the system to achieve the necessary 
resonant frequency. 


4) The describing function technique may be readily 
employed as an analytic tool in the study of a 
missile autopilot system having a relay servo. 


12 IRE TRANSACTIONS ON AUTOMATIC CONTROL 


APPENDIX 
DERIVATION OF THE DESCRIBING FUNCTION 


Considering the waveform to be a periodic odd func- 
tion, with the period divided equally, the time function 
in terms of a Fourier series is 


f® = DX q sin n(wt — a) (9) 


where 
ee (oe , 2F 
an = =f f(x) sin nxdx = —(1—cosmm). (10) 
Tv 0 TH 


From (9) we have that one term of f(t) is: 
A, cos nwt + B, sin nwt 
where 
An = — Gp» SiN na. 


By = Gn COS Na. 


Consider initially the effects of harmonics that are of 
a higher order than the fundamental as being negligible; 
4.e.,n=1 only. 

Substituting into (10) we find that 


4F 
C= bi 
T 
and 
BV ae 
An sin a 
Tv 
4F 
By, = — cosa. 
Tv 
However, from Fig. 5, 
, h 
sin a = — 
x 
and 
4Fh 
1 = 
1x 
4F~/x? — }? 
1x 
The describing function is 
VAP Be 


N(x) = (11) 


November 


and 


—1 x 


(12) 


ie ats ha [ep 
Nie) n/Av + Be 


where 
B = ea Bi/ Ai. 


Upon substituting for A, and B,, we arrive at the 
required relationship that 


(13) 


= a ff 8 : 
= =(<), /—180° + sin~! h/x. 
NG wee 


GLOSSARY 


Roman Symbols 


A =roll rate gain (pounds/rad/sec). 

A,, An=Fourier coefficients. 

@, @ = Fourier coefficients. 
B=roll attitude gain (pounds/rad). 
F, f=control force (pounds). 
h=hysteresis (pounds). 
Iz, Iy, Iz, Ire =moments and product of inertia (slug 
—feet?). 
j =jet control subscript, also ~/—1. 
L=roll moment (pound — feet). 
N(x) =servodescribing function. 
n =integer of harmonic. 
P=roll rate (rad/sec). 
Q=pitch rate (rad/sec). 
R=yaw rate (rad/sec). 
s =Laplace operator. 
t=real time (sec). 
X =amplitude of servoinput signal, roll axis 
of missile. 

Y=missile pitch axis. 
Z=missile yaw axis. 


Greek Symbols 


a= phase lag. 

6 =phase angle (deg). 

T =system time delay (sec). 

@=roll angle (rad). 

w =radian frequency (rad/sec) 

¢ =damping factor (nondimensional). 
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Root Locus Method of Pulse Transfer Function for 


Sampled-Data Control Systems” 
MASAHIRO MORIt+ 


Summary—The pulse transfer function is a powerful tool for 
analysis and synthesis of sampled-data systems. The author has ex- 
tended the root-locus method! to the region of the pulse transfer func- 
tion so that we can analyze and synthesize the sampled-data systems 
by the analogous method for the continuous-data systems. 

The results are as follows: the rules by which the root loci of the 
pulse transfer function can be plotted are identical to those of the 
continuous-data systems. The conditions for absolute and relative 
stabilities for the sampled-data systems on the z plane are conformal 
to those for the continuous-data systems on the s plane. 

Existence of a finite settling time in a sampled-data control sys- 
tem has been pointed out by the root-locus method. 


INTRODUCTION 

( ENERALLY, the pulse transfer function (ptf)?? 

G* is a rational function of e*? and can be written 
N*(e*7) N*(z) 

D*(e*7) ee D*(z) - 


ater G*(z) (1) 
where z=e*?, T is the sampling period, and N*, D* 
denote polynomials in zg. The order of z in the denom- 
inator D* is equal or higher than that of the numerator 
N*. G* isa function of s, hence, of z. Therefore, the root 
locus may be defined on the s plane as well as on the z 
plane. But, since G* is a rational function of z and a 
transcendental function of s, it is more expedient to de- 
fine root loci in the z plane than on the s plane. 

The root loci of pulse transfer functions are defined as 
plots of the variations of the poles of the closed-loop 
pulse transfer functions on the z plane with changes in 
the open-loop gain. In this paper, the root loci are 
termed “z loci,” and correspondingly, the root loci on 
the s plane for the continuous-data systems are termed 
ag slog.” 

With reference to Fig. 1, the closed-loop ptf G,.*(z) is 
given as 


C*(@) * KG*(z) 
R*(z) 1+ KG*(z) 


where G*(z) is the open-loop ptf, that is, G*(z) is the z 
transform of impulsive response of the system whose 


G*(z) = (2) 


* Manuscript received by the PGAC, January 1, 1956. Presented 
at the Symposium of The Society of Automatic Control in Japan, 
June 29, 1954. Published in Automatic Control (Japan), vol. 2, pp. 
5-12; March 1, 1955. 

+ Univ. of Tokyo, Yayoiché, Chiba City, Japan. 

1W. R. Evans, “Graphical analysis of control systems,” Trans. 
AIEE, vol 67, pp. 547-551; 1948. P 

2 F, Hurewicz, “Theory of Servomechanisms, H. M. James, N. B. 
Nichols, and R. S. Phillips,” M.I.T. Rad. Lab. Ser., McGraw-Hill 
Book Co., Inc., New York, N. Y. vol. 25, ch. 5; 1947. ; ged 

3 R. H. Barker, “The pulse transfer function and its application 
to sampling servo systems,” Proc. IEE, part 4, monograph no. 43, 
pp. 302-317; July, 1952. 
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Fig. 1—Typical sampled-data control system. 


transfer function is Ga(s)Gr(s)G,(s). K is the open-loop 
gain. 

Accordingly, the z locus for the system shown in Fig. 
1 is the locus of the root of the following characteristic 
equation as K is varied from zero to infinity. 


1+ -KG*(z) = 0. (3) 


RULES FOR CONSTRUCTION OF 2 LOCI 


As mentioned above, the continuous system transfer 
function is the rational function of s, and the ptf is the 
rational function of z. Accordingly, the rules of the s 
loci which are based upon the rational form must be 
applicable to the z loci. Such rules are as follows: 

1) For a gain K=0 the roots of (3) coincide with the 
poles of G*(z), and for a gain K= the roots coincide 
with the zeros of G*(z). Thus, the z loci start from the 
poles and arrive at the zeros. 

2) The z loci are always symmetrical about the real 
axis. 

3) Phase condition: From (3), 


KG*(g) =, 1. 


Hence, with reference to Fig. 2, the phase condition is 
given as 


> B— do a= — 180° + 360° X x (4) 


where @ is the angle of the vector from a zero of G*(z) to 
the root Q, and a is the angle from a pole of G*(z) to Q 
as shown in Fig. 2, and 7 is an integer. 

4) Asymptotes: If the number of poles m exceeds the 
number of zeros m, there are n—m branches which 
terminate at the zeros located at infinity. The directions 
of the infinite asymptotes of the zg loci are shown in 
Table I. All asymptotes intersect at a point pa on the 
real axis, where 


n m™m 
DS Pi — D8 
4=1 i=1 
Pa 
nw—-m 


14 IRE TRANSACTIONS ON AUTOMATIC CONTROL 


10 Real 
Pi 
X: Pole 
O: Zero 


Fig. 2—Illustration for gain and phase conditions. 


TABLE I 
DIRECTIONS OF ASYMPTOTES 


n-™ Directions 
1 e 180° 
2 —90° +90° 
3 —60° +60° +180° 
4 —45° +45° —135° Snes 


provided p;=pole and z;=zero. 

5) Angles of departure and arrival: The angles at 
which the z loci leave the poles and arrive at the zeros 
can be evaluated from (4). For example, in Fig. 3 we 
wish to determine the angle departure of the z locus 
from pole P. Let us consider a point Q on the z locus a 
small distance away from the pole P. This point Q must 
be at an angle a; from P such that, 


(8: + Bo) — (a1 + a2 + a3) = — 180° + 360° X xn. 


6) The z locus exists on a part of the real axis where 
an odd number of poles plus zeros are found to the 
right. This may be verified easily by (3). 

7) Intersection with the real axis: If the displacement 
from the real axis is A and this is intentionally made 
very small, then in Fig. 4, 


180° — a, = A/ 


az = A/Iy 
bi = Ay ls. 
Applying the condition to rule 3), 
Bi — a2 — a, = — 180° 
MA eA): 
Thus, 
ie ie ctr h: (5) 


Thus the intersection may be determined as a point 
where the equality indicated in (5) exists. 
8) Gain condition: From (3), 
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Fig. 3—Evaluation of angle departure. 
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Fig. 4—Determination of intersection with the real axis. 


1 
Ste 
| G*(e) | 
A graphical procedure is convenient to evaluate values 


of K along the z loci by this relation. For example, if 
G*(z) =0.6z/(z—1)(z—0.4), with reference to Fig. 2, 


EL VE=04)| _ POX PAO 
0.62 0.6 X 00 
1.34 X 0.84 
ca IN 0.630 


Thus, gain K can be evaluated in terms of the vectors 
from the poles and zeros to the point on the zg locus. 

9) Centroid of the roots: As mentioned by Yeh,‘ the 
centroid p of the roots, which is given by 


1 n 
[oy > Pi 
WM ¢=1 
is an invariant of the system if »—m=2; it will shift by 
an amount equal to K/n if n—m=1. 
10) Loci of simple ptf: The ptf having m poles and m 
zeros shall be referred to as type (n, m). Loci of type 


(22, mS2) are straight lines or circles as shown in 
Table II.4 


*V. C. M. Yeh, “The study of transients in linear feedback sys- 
tem by conformal mapping and root-locus method,” Trans. ASME, 
vol. 76, pp. 349-361; 1954. 
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TABLE II 
Z LOCI OF SIMPLE PULSE TRANSFER FUNCTIONS 
Type Shape 
(1, 0) Straight line from the pole to minus infinity on real axis 
leet) Straight line between the pole and the zero on real axis For type (2, 1) and type (2, 2), 
i, pole=i+ Jv, u2—jv 
For two real roots: Segment between the two poles and Zero=&-+ jn, &2—jn 
(2, 0) perpendicular bisector of the segment Radius =r 
For conjugate complex roots: Straight line through the 
two roots (except segment between the roots) 
Qxi) Circle (its center is on real axis) Center = & 7? = (1 — €1) (u2— £2) +r? 
(mimetv?) — (£1& +7?) 
Center, 0= 
; (m1 toe) — (£1-+ £2) 
r2) Circle Sag he (uaF Me) (Mamet v?) — (&1+ £2) (Era-+7?) 
= 
(uit M2) — (E1-+ £2) 


STABILITY CRITERION 


Generally speaking, if all singular points of the ptf 
are located inside the unit circle, the sampled-data sys- 
tem is said to be stable with regard to output values 
only at sampling instants. If at least one singular point 
lies outside the unit circle, the system is unstable. 

The system is stable if all roots of (3) are located in- 
side the circle. This constitutes the necessary and suffi- 
cient condition of the closed-loop system stability. 

The unit circle in the zg plane corresponds to the 
imaginary axis of the s plane. 


TRANSIENT RESPONSE 
Relation Between the Location of the Root and the Shape 
of the Indicial Response 


Information on the effects of the roots on the transient 
responses are given by inverse z transformations. When 
the reference input is the unit step, the output sequence 


where Z—! denotes the inverse z transformation.’ Thus, 
in general, 


(6) 


ae zfore 


d f 
KIT @—+) It@- &)? +22} 


1 
Cy = 5 ; h 
21 iat Cea II! (g — wi)? + v2} 
_* _ ge-adg (7) 
(g =~ 4} 


where the path of integration in (7) is a closed contour 
enclosing all the singular points of G,*(z) and r;=real 
root of (3), ustjv;=conjugate complex roots of (3), 


5 C. H. Smith, D. F. Lawden, and A. E. Bailey, “Characteristics 
of Sampling Servo Systems,” in A. Tustin, “Automatic and Manual 
Control,” Butterworths Scientific Publications, London, Eng.; 1952. 
j ee. 


z;=real zero of (2), £;+jyni=conjugate complex zeros 
of (2), and g+2h2d+2f. 
Evaluating the contour integral, we obtain 


De | R: 
saan 'k 


rj 


(ui? + v:7)*/? cos (RO; + $3) (8) 
where R; is residue of the real root 7;, I; is residue of 
the complex root u:tjv:, 0; is phase of the complex root 
Mityjvi, and ¢; is phase of the residue /;. 

In the following discussions, it is assumed that all 
roots have absolute values less than unity. In (8), if 7; 
is positive, r;* decreases uniformly as k increases, and 
converges to zero. If 7; is negative, 7;* is oscillatory and 
converges to zero. Hence, a positive real root represents 
a uniform variation and a negative real root represents 
a damped oscillation. 

Now, let us focus attention on the second term. It is 
apparent from previous assumptions that the abso- 
lute value of a pair of conjugate complex roots, 7.e., 
(u;?+y,7)1/? is less than unity. Therefore, each com- 
ponent of the second term which comes from a pair of 
conjugate complex roots describes a damped oscillation. 

Consequently, it can be concluded that a uniform 
variation comes from a positive real root and a damped 
oscillation is from a negative real root or from a pair of 
conjugate complex roots. 

Also, it is clear from above that the intersection of 
the z locus with the positive real axis indicates critical 
damping. 

In addition to these results, it is noted that output at 
the steady state comes from a real root at the point of 
(+1, 0). This can be deduced from the first term of (8) 
or from the final value theorem for z transforms.® 

These relations are conformal to those at the continu- 
ous-data systems in the s plane, and are arranged in 
Table III. 


6 J. R. Ragazzini and L. A. Zadeh, “The analysis of sampled-data 
systems,” Trans. AIEE, vol. 71, part 2, pp. 225-232; 1952. 
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TABLE III 
CONFORMAL RELATIONSHIP BETWEEN THE s PLANE AND THE 2 PLANE 


s plane 


z plane 


Origin 


(+1, 7:0) 


Real axis 


Positive real axis 


Imaginary axis 


Unit circle 


Right half plane (width is 27/T) 


; Outside of the unit circle 


Left half plane (width is 27/T) 


Inside of the unit circle 


Stability criterion Left half plane 


Inside of the unit circle 


One uniform variation One negative real root 


One positive real root 


One damped oscillation | One pair of conjugate complex roots 


One negative real root or one pair of conjugate complex roots 


Critical damping 


Intersection of s locus with negative real axis 


Intersection of z locus with positive real axis 


Output at steady state | Root at the origin 


Root at point (+1, 7-0) 


Constant —f curve Straight lines /, 1’ in Fig. 5 (a) 


Heart-shaped curve in Fig. 5 (b) 


Constant-¢ Curve 


Now the damped oscillation discussed above will be 
studied more quantitatively in the following. 

An output sequence from a component of the second 
term of (8) is, 


c, = 2| 1; 


| g:|* cos (0: + 43), (9) 
where g;=itjvi, and | gi = (p2+v,?)1/?. 

If one assumes a continuous output wave of damped 
cosine form, c, being the sampled values, one can apply 
the concept of damping, and define the following am- 
plitude ratio: 


im (Gp) hee dn 10; ; (10) 


(Ck) k=0 


The relation between rz and the complex root q; is ob- 
tained by substituting (9) in (10) as follows: 


oy, 
Vd 


py eea (38 6” 
2| I:| | 


| 9s 20". 


(11) 


° cos g; 


Hence, the condition for an output with amplitude 


ratio rq is: 
( gi 
bs gi 


The locus of (12) is a logarithmic spiral as is shown in 
Table IV and Fig. 5(b). 

On the other hand, the same damping condition for the 
continuous-data system on the s plane is a line / in Fig. 
5(a). Its equation is 


raiil em 


(12) 


6; 


e2tale =F 4; 


(13) 


where s =o-+]w. 


TABLE IV 

ABSOLUTE VALUE |g;| AND PHASE 6; FOR ¢=0.2 (r=0.25) 
Phase 6; Absolute value Phase 6: Absolute value 

Ig! lg| 

c 1.00 51/8 0.65 

hs 0.92 31/4 0.59 

ae 0.84 7/8 0.55 

3n/8 0.77 x 0.50 

a/2 0.71 
je"sinw 


Const- line 


(a) (b) 


Fig. 5—Conformal relationship between the s plane and z plane. 


Therefore, 


(14) 


a/w = — log ra. 
2a 


Transforming this line on the z plane, we get 
est i er l eieTl = fg ek ieee 
which has the following complex value: 


Gem value = e°F = 7yeTl2n 


wT, (15) 


phase 
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Comparing (12) with (15), one can see that the spiral in 
Fig. 5(b) is conformal to the straight line ] in Fig. 
5(a). Thus, the constant amplitude ratio curve for the 
sampled-data systems coincides with the transform of 
that for the continuous-data systems. Also, it is deduced 
that w, the imaginary part of s, is angular frequency of 
the output cosine wave. 

The above discussions based upon rg can be modified 
in the following statements in terms of conventional 
damping ratio ¢. Using the angle 8 shown in Fig. 5(a), 
¢ can be written as: 


(== Sin (oh. (16) 
Referring to (14), we have 
o 
tan B = — = — log rq. (17) 
e) Qr 
Thus, 
: 1 
¢ = sin {tan (— log rt ’ (18) 
2a 
Therefore, (12) can be rewritten as 
| gi] = eG/v Ie 0: 
( (19) 
Z4i = 6;. 


Though G,*(e:7) contains poles of infinite number, 
periodic on the s plane, this is equivalent to G,*(z) 
which is characterized by a finite number of poles on 
the z plane. By the conformal relationship of s and z 
(see Fig. 5), the inversion integral (7) is specified 
uniquely by choosing a contour about one periodic strip 
on the s plane. Consequently, it is obviously justifiable 
to choose this contour about the first strip bounded by 
the lines +jr/T. Therefore, though the logarithmic 
spiral continues infinitely towards the origin as 6; in- 
creases, it should be limited to the part of 0;S7. 

Thus, the heart-shaped curve in Fig. 5(b) is the con- 
dition for constant ¢. If the conjugate roots are located 
inside the heart-shaped curve the damping is stronger 
than that of this condition, and if the roots are located 
outside, the damping is weaker. Table IV shows polar 
coordinate of the constant-¢ curve for ¢=0.2 (or 
ra=1/4). 

The heart-shaped constant-f curve has an intersection 
with negative real axis for 6;=a. Thus there exists a 
negative real root which gives a damped oscillation of 
the specified ¢ value. This fact is not surprising because 
a negative real root inside the unit circle causes a 
damped oscillation as mentioned before. 


RELATIONS BETWEEN z Locus AND s Locus 


It should be noted that the root loci on the gz plane, 
of sampled-data systems, 1+G*(z) =0, are not the con- 
formal transformations of root loci of 1+G(s) =0 on the 
_ s plane. 
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As pointed out by Ragazzini and Zadeh,‘ the follow- 
ing relation was discovered by Poisson. 


1 


i > G(s + 2n1j/T), 


—— 0 


Gals (20) 


where 7 is an integer. 
Therefore, the z locus of (3) is conformal to root locus 
on the s plane given by the following: 


ce} 


1 
1+ = >> G(s + 2nrj/T) = 0. 


n>=—o 


(21) 


However, plotting this root locus on the s plane is diffi- 
cult. 


DEAD TIME 


Handling a dead time, LZ, by the s locus of a continu- 
ous-data system’ is not simple, because its transfer func- 
tion has term e~*” in the numerator. 

On the contrary, the z locus of a sampled-data system 
with a dead time can be plotted as easily as the system 
with no dead time (see Fig. 6), because the ptf Gz*(z) 
which has a dead time L is a rational function of z as 
follows: 

1) When L=xT, 

Gi*(z) = G*(z)/2" (22) 
where G*(z) is ptf of the system having no dead 
time, and 7 is positive integer. 

2) When L=S6T, 


Gi*(2) = 27 D0 gerimre* 


k=0 


(23) 


where gir is weighting sequence of the system 
whose ptf is G*(z), and 0S$6<1 and m=1-—6. 


FINITE SETTLING TIME 


It is evident that in the continuous-data control any 
linear system has an infinite settling time (time to reach 
zero error after applying a step input). 

On the contrary in the sampled-data control, when 
the z locus passes through the origin [see the right sides 
of Fig. 6(a) and (b) |, itis possible to cause the system to 
have a finite settling time® so far as the controlled vari- 
able at the sampling instants is concerned. If the gain is 
equal to the value of K at the origin, the impulsive re- 
sponse derived from the zero root becomes a finite poly- 
nomial in z~!. This will be demonstrated in example 2). 


7A. Nomoto, “Contribution of the root-locus analysis of the 
feedback control system,” Proc. Second Japan Natl. Congress for 
Appl. Mech., pp. 359-362; 1952. 

8 R. C. Oldenbourg, “Deviation Dependent Step-by-Step Control 
as Means to Achieve Optimum Control for Plant with Large Dis- 
tance-Velocity Lag,” in A. Tustin, “Automatic and Manual Control,” 
Butterworths Scientific Publications, London, Eng.; 1952. 
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4 Imaginary * Imaginary x Mee Ck 
Reference \ bait 
trolle 
z-Plane nitro! 
Zero Order Controlled Variable 
Hold System 
ef ! Real 
Fig. 7—Sampled-data control system used in illustrative 
example 1). 
™ 
1+K Sop ) 
AT maginary 
z-Plane 
dcosé8 
cos ™é 
amen L= 10 sec. 
Fig. 8—Indicial response of the controlled system used in 
illustrative example 1). 
the loci start from the origin and (+0.78, 0) and one of 
Vane z*-2dcosP _ fin gr Zcosmb = dcosdP the loci terminates at (—1.75, 0). From rules 2) and 4), 
© 2'-22d cosp +d° z'-2deos@ + d another locus terminates at minus infinity. From rule 6), 
(b) the parts of the real axis between the origin and (+0.78, 
(a) Gs) = (b) G(s) = Se 0), and between minus infinity and (—1.75, 0) con- 
(s + a) (+a)? +8 stitute sections of the loci. With reference to Table II, 
Left ee without dead time, characteristic equation is: since the type of this ptf is (2, 1), the other part of the 
1+KG*(z)=0, 


Right sides: with dead time less than 7, characteristic equation 
is: 1+KG*z(z) =0. 
[See¥(23).] 


Fig. 6—Examples of z locus. 


Illustrative Examples 


Example 1): In the first example, the gain adjustment 
of the system shown in Fig. 7 is considered. It is as- 
sumed that the specification for damping is ¢=0.2, the 
continuous transfer function of the controlled system is 
e-*4/(sr+1) and the time constant 7 is 60 seconds, the 
dead time L is 10 seconds. It is postulated also that the 
sampling period T is 15 seconds and the order of the hold 
circuit is zero, that is, its transfer function is (1—e7*7)/s. 
(See Fig. 8.) 

The characteristic equation is written as 


1+ Kf ae it ah 24 
z (g — e-T!7)z Ties: “> 
Hence, 
0.08(z + 1.75 
i+ K COB ah {9)) = (25) 


2(z — 0.78) 


The z loci can be drawn on the basis of the rules from 1) 
to 10) as follows: the zero is (—1.75, 0), and the poles 
are the origin and (+0.78, 0).. Therefore, from rule 1) 


loci is a circle with a radius 2.11, and the center at 
(—1.75, 0). Thus the z locus can be drawn easily as 
shown in Fig. 9. 

Constant-f¢ curve for ¢=0.2 can be obtained by plot- 
ting the complex value q; of (19). 

Then the gain K for ¢=0.2 is determined by finding 
the value at the intersection of the z locus with the con- 
stant-¢ curve. Absolute values of three vectors are re- 
quired for the evaluation, and by rule 8). 


_ | 2(e — 0.78) 0.76 X 0.90 
ST 0)08(¢.5. 1:75)! 0.08 X 247 


Under this condition, the indicial response for the 
reference input may be obtained by (7). [For this pur- 
pose, it is more convenient to expand the integrand of 
(7) into a power series in z~!. The coefficient of the term 
of z~* is equal to the value of the output at the kth sam- 
pling instant c,].° 


Che 2X O05 oe. 0e* 


X cos {(k — 1)71.5° + 45°} +.0.80. (26) 


Fig. 10 shows this response. 
Example 2): Here a compensation of a sampled-data 
control system is demonstrated. The ptf of the system 


* J. R. Ragazzini and A. R. Bergen, “A mathematical technique 
for the analysis of linear systems,” Proc. IRE, vol. 42, pp. 1645- 
1651; November, 1954. 
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Fig. 9—Z locus of the sampled-data control systems used in 
illustrative example 1). 


K T=0 
Real 


(b) 


Fig. 11—(a) Sampled-data control system used in illustrative 
example 2), without compensation. (b) Its 2 locus 


—-——(J-z"') aa 
Delay line 
En-I 
Out put 


en " (nee) 


Fig. 13—Compensating pulsed network used in illustrative exam- 
ple 2); en is sampled value of input waveform e(t) at time nT. 


& : Double pole 
O'"Zer0 


Reference input 


O5+ 


Gmetieet3t 4hest ci 


Fig. 10—Indicial response of the sampled-data control system 
used in illustrative example 1), when gain, K, is 4.0. 


Controlled 
system 


(b) 


Fig. 12—(a) Sampled-data control system used in illustrative 
example 2), with compensating network. (b) Its z locus. 


— 


n 


Contyolled 
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— 


\ sot ing time 


Fig. 14—Indicial response for unit step reference input of the 
system shown in Fig. 11, (a) KT=1.0. 
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which has the continuous transfer function 1/s? is, 


G*(z) = 2T/(z — 1)? (27) 
where T is the sampling period. 

If we construct a sampled-data feedback system sim- 
ply as shown in Fig. 11 (a), without any compensation"? 
the characteristic equation will be, 


foateste 0 (28 | 

= CET ) 
Fig. 11(b) shows this z locus. According to the stability 
criterion, this system is not stable for any value of KT. 

The system may be stabilized by adding a compensat- 
ing unit in series with the controlled system. The com- 
pensating unit is a pulsed network as shown in Fig. 13, 
and its output at time 77 is the difference of two sample 
values of input at mT and (~—1)T7,'! where x is an in- 
teger. Then the ptf of this network is 1—z7!. Thus, the 
characteristic equation becomes 


£ 


Bs 


Fig. 12 shows the compensated feedback system and 
its z locus. The system becomes stable for 0<AT <2. 

If K has the value at the origin, the system has the 
finite settling time. That is, for this value of K the 
closed-loop ptf is, 


10 In this case, the input to the controlled system is an impulse 
since the system has no hold circuit. 

u W. K. Linvill and J. M. Salzer, “Analysis of control systems 
involving digital computers,” Proc. IRE, vol. 41, pp. 901-906; 
July, 1953. 
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(eel) 
Ge 

(a> 1) 


Fig. 14 shows the control action under this condition. 
The settling time is only one sampling period and the 
system has no overshoot.” 


= 


G.*(z) = (30) 


Co 


KT=1 


DISCUSSION 


In short, the rules for graphical construction of z loci 
are equivalent to those for s loci. 

Absolute and relative stabilities can be deduced from 
conformal mapping conditions between s and zg. The 
s plane may be transformed to the z plane in many folds. 
However, the pertaining region of the s plane should be 
limited to the strip +j7/T. 

When the z locus passes through the origin, it is pos- 
sible to cause the sampled-data system to have a finite 
settling time. 

The author would like to emphasize the usefulness of 
the root-locus method not only in- continuous-data 
systems but also in sampled-data systems. 
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The z transform method has the limitation of yielding informa- 
tion only at the sampling instants. Therefore, the concept of the 
finite settling time should be limited only to the controlled variable 
at the sampling instants. But, as a matter of fact, the system in this 
example has a finite settling time T in a continuous sense. That is, 
after the time 7, the error of the system becomes zero not only at 
sonal oa a: instants but also everywhere as shown by dotted line 
in Fig. 14. 


GS eo 


1957 


IRE TRANSACTIONS ON AUTOMATIC CONTROL 21 


Input-Output Analysis of Multirate Feedback Systems’ 


GEORGE M. KRANCT 


Summary—A general analytical technique described in this 
paper permits the extension of Z-transform methods to sampled- 
data systems containing synchronized switches which do not operate 
with the same sampling rate. Sampling periods of each switch are 
first expressed in the form 7/p; - - - T/p, (where p, - - - p, are inte- 
gers not equal to zero) and then it is shown that each switch with a 
period 7/p can be replaced by a system of switches and advance and 
delay elements where each switch operates with a sampling period T. 
In this way, the original sampled-data system can be represented by 
an equivalent system containing switches operating with the same 
sampling rate. The general solution of such equivalent systems is out- 
lined in this paper. 


INTRODUCTION 


N THE modern technical literature, a great deal of 
l attention has been paid to the analysis and syn- 

thesis of feedback systems operating with discrete 
or sampled data.!~® The great majority of the systems 
described in the literature have one feature in common: 
they contain a number of switches which operate in 
synchronism with the same sampling rate. Such sys- 
tems will be referred to as single-rate sampled systems 
to be distinguished from multirate sampled systems, 
defined here as sampled-data systems containing syn- 
chronized switches which, however, do not operate with 
the same sampling rate. 

The practical application of multirate systems arises, 
in many cases, out of problems connected with the anal- 
ysis and synthesis of standard, single-rate sampled- 
data systems. Several investigators have used such sys- 
tems as convenient mathematical analogs. The results 
they obtained were largely by-products of the work done 


* Manuscript received by the PGAC, December 4, 1956. This 
research was supported by the U. S. Air Force under Contract No. 
AF 18(600)-677. Presented at WESCON, Los Angeles, Calif., August, 
1956. 

» + Dept. of Elec. Eng., Columbia University, New York 27, N. Y. 
F~ 11. A. MacColl, “Fundamental Theory of Servomechanisms,” 
D. Van Nostrand Co., Inc., New York, N. Y.; 1945. 

2W. K. Linvill, “Sampled-data control systems studied through 
comparison of sampling with amplitude modulation,” Trans. AIEE, 
vol. 70, pp. 1779-1788; February, 1951. : 

3J, R. Ragazzini and L. A. Zadeh, “Analysis of sampled-data 
systems,” Trans. AIEE, vol. 71, pp. 225-234; November, 1952. | 

4A.R. Bergen and J. R. Ragazzini, “Sampled-data processing 
techniques for feedback control systems,” Trans. AIEE, vol. 1s) 
pp. 236-247; November, 1954. ‘ : By 

5 R. H. Barker, “The pulse transfer function and its application 
to sampling servo systems,” Proc. IEE, vol. 99, pp. 302-317; Decem- 

Was 
ee o K. Linvill and R. W. Sittler, “Extension of conventional 
techniques to the design of sampled-data systems,” 1953 IRE Con- 
VENTION RECORD, part 1, pp. 99-104. : 

7G. V. Lago and J. G. Truxal, “The design of sampled-data 

feedback systems,” Trans. AIEE, vol. 73, pp. 247-253; November, 
54, 

o 8 J, Sklansky, “Network Compensation of Error-Sampled Feed- 

back Systems,” Ph.D. Dissertation, Dept. of Elec. Eng., Columbia 

University, New York, N. Y.; 1955. ; 

9G. M. Kranc, “Compensation of an error-sampled system by a 
multirate controller,” Trans. AIEE, vol. 77, pp. 149-159; July, 
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in other fields of sampled-data theory. Thus, in order to 
evaluate the output between sampling instants, Linvill 
and Sittler® introduced a fictitious switch operating at 
twice the sampling rate of the system. Essentially the 
same approach was later adopted by Lago and Truxal.’ 
Sklansky® has shown how an error-sampled feedback 
system can be compensated using continuous networks, 
by first constructing an approximate model which con- 
tained an additional switch operating at a rate equal ton 
times the sampling rate of the system, ” being an in- 
teger. He obtained an expression for the closed-form 
Z transform of the output of the approximate system. 
The method he used extends ideas of his predecessors 
to the case where two switches, operating at sampling 
rates related by an integer, are included in the feedback 
loop. 

In addition to the cases where multirate sampling is 
used for purely computational purposes, there is a 
possibility of the actual physical existence of such sys- 
tems. This may occur in the case of flight-control sys- 
tems, when a data link computer operates at a rate dif- 
ferent from the speed of the radar antenna; or in a situa- 
tion where multirate sampling is used purposefully to 
improve the performance of a given system. A sampled- 
data processor or digital controller operating at a higher 
rate than the basic sampling rate of the system would 
exemplify the latter case.° 

Although sampled-data systems belong to the class 
of time varying systems, it is possible to treat them very 
effectively from the point of view of time-invariant or 
fixed systems through the use of Z transforms. It is 
shown in this dissertation that multirate sampled sys- 
tems can be treated also from the point of view of fixed 
systems, using standard tools of sampled-data mathe- 
matics. In this way, the standard methods of Z trans- 
forms are extended to cover the case of multirate 
systems. 


BRIEF SURVEY OF THE FIELD 


_It has been mentioned already that some theoretical 
investigations of multirate systems have been carried 
out in connection with problems arising out of the an- 
alysis and synthesis of single-rate systems. The results 
of these investigations presently will be derived again 
in a manner which is consistent with the general ap- 
proach of this paper. In this way the notation, which 
is largely due to Sklansky, and the basic problem will 
be introduced. In all examples in this and other sec- 
tions, switches are assumed to first close simultaneously 
at time t=0. 

Considering first the system. shown in Fig. i(a), 
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Fig. 1—Multirate open loop systems. 


where the switch S; operates with the sampling period 
T, it can be written that 


C(s) = H(s) E(z) (1) 


where C(s) is the Laplace transform of the continuous 
output c(t), H(s) is the Laplace transform of the im- 
pulsive response of the network H, E(z) is the Z trans- 
form of the input e(t), and z=e?. The Z transform of 
e(t) can be denoted also by Z[e(t)| or Z[E(s)] where 
E(s) is the Laplace transform of e(#). Since E(zg) is 
periodic with respect to w)=27/T, Z transforming (1) 
produces the well-known result: 


C(z) = Z2[C(s)] = Z[HG)E@] = HE) (2) 


where H(z) =Z|H(s)] is referred to sometimes as the 
“pulsed transfer function.” C(z) represents the output 
sequence of the system shown in Fig. 1(a) where the 
switch S, operates with the same sampling... period 
... / as the switch Sj. 

Let it be assumed now that the period of the switch 
Sis reduced from 7 to T/n, n being an integer. This is 
shown in Fig. 1(b) where the output is sampled at the 
rate n times higher than the sampling rate at the input. 
The corresponding Z transform will be denoted by Z, 
and the output sequence is given by Z,,[C(s) |. Z, trans- 
formation of (1) gives 


Z,{C(s)] = Z,[H(s)E@)]. (3) 


Let C(zn) 4 Z,[C(s)] and H(z,) 4 Z,,[H(s) ]. Then, since 
E(g) is periodic in wo, and must also be periodic in nwo, it 
can be written that 


C (8n) = H(én)E(z) = H(2n)E(2n") (4) 


where 2, is defined as z!/"=e°7/”, It is important to note 
that C(z,) is not equal to C(z) with z replaced by z,, and 
that since H(z,) represents the impulsive response of H 
sampled with the period T’/n, it can readily be found 
from the tables of Z transforms simply by replacing all 
T by T/n, and all z by z, in the expression for H(z). On 
the other hand, E(z,”) is obtained by replacing every zg 
in the expression for E(z) by 2,”. 
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In the system shown in Fig. 1(c) the switches S; and 
S2 both operate with the same period T’/n. Therefore, it 
follows from (1) that 


C(an) = H(2n)E(Zn) (5) 


where E(z,,) 4 Z,[E(s) |] is obtained in the way previously 
indicated. To obtain the output transform C(g), it is nec- 
essary to note that only every mth sample of the se- 
quence C(z,) can be transmitted through the switch Sz, 
which operates with the sampling period 7. Let 
P,,[C(zn) | 4 every nth sample of C(zn). Then 


Cz) = Z[C(zn)] = PalC(en)] = Pal E(@n)H(En)]. (6) 


Sklansky? has shown that, in general 


P,[C(En)] = oP y Clan Orkl™, (7) 


WN k=1 


and therefore, 


done : 
Ce = = »>y E(Snei 27F!™)) A (Zyeh 2tkIn)) | (8) 


W K=1 


It is interesting to note that in (8) the output trans- 
form C(z) does not appear any more as a product of the 
input transform and the pulsed transfer function. In 
fact, the concept of the pulsed transfer function has no 
meaning in this case unless a time dependence is intro- 
duced. Since the switch S: operates with the rate equal 
to m times the rate of the output samples, its presence 
is entirely superfluous and therefore it can be replaced 
by a short circuit without any effect on the output 
transform C(g) as shown in Fig. 1(c). Thus, one can 
treat the diagram of Fig. 1(c) as representing a conven- 
ient mathematical model of the system in Fig. 1(d). 

The procedure outlined above now will be applied to 
the feedback system shown in Fig. 2(a). The diagram is 
redrawn in Fig. 2(b) with a superfluous switch, operat- 
ing with the period 7T/n, inserted after the network H. 
It follows from this diagram that 


E(z) = R(z) = Z| H(en)C(Zn)]. (9) 


On the basis of results previously developed for the 
open loop case it can also be stated that 


ClEn) = Gl2n) E(2n”). (10) 


Next, eliminating C(z,) from (9) and (10) and making 
use of the periodic properties of E(z,”")=E(z), it can be 
shown that 

R 
ce (z) 
1+ Z[H(z,)G(z,) | 


Since Z| H(2n)G(zn)] implies P,, [H(zn)G(zn) ], one can 
make use of the identity stated in (7); also, E(z) and 


R(z) can be replaced by E(z,”) and R(zn"), respectively. 
Eq. (11) can, therefore, be restated as 


(11) 
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E (z) 


C(s) 


(s) + 


Cs) 


C (Z,) 
oy UI a 


(b) 
Fig. 2—Multirate feedback system; sampling rates related by 
an integer. 
R22" 
E(@,”) = (12) 
= ; 
1 ae H (zn€227*!") G(gyet27*/”) 
NM k=1 
Eliminating E(z,”) from (10) and (12), one obtains 
R(En")G(En) 
8 ere (13) 
fet > H (2,227!) G(z,¢827#/7) 
N k=1 
Let 
AaC(Z 5) 
Kz) = ) 
R(Z,.”) 
then 
G(Zn) 
K(2n) = (14) 


n 


tee > H (Zne?27*!") G (Zne227*/”) 
N k=1 

K(zn) can be interpreted here as the Z, transform of the 
impulsive response of the feedback system shown in 
Fig. 2(a). It follows from (14) that K(z,) can be ex- 
pressed as a ratio of two polynomials in z,~! and there- 
fore in order to determine whether the system shown in 
Fig. 2 is stable it is sufficient to examine the roots of the 
denominator. In analogy to the single-rate sampled 
systems, K(z,) is a stable transfer function if its denom- 
inator vanishes only for values of z, corresponding to 
points within the unit circle on the z, plane. 


INTRODUCTION TO THE SWITCH 
DECOMPOSITION METHOD 


The procedure outlined in the previous section is 
limited in its application to cases where a sampled-data 
system contains switches operating at sampling rates 
‘whose ratios can be expressed by an integer. A different 
analytical approach will now be outlined, whose main 

advantage is that it does not suffer from the above men- 
tioned limitation. Let F(s) be the Laplace transform of 
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Fig. 3—Switch decomposition diagram. 


the continuous input f(t) to a switch, operating with 
the sampling period T/n. Then the output sequence is 
given by Z, transform of f(t), namely 

F(2n) = >) f(RT/n)en*,  k=0,1,2,- (15) 

k=0 

The same output sequence is obtained from an analo- 
gous mathematical model consisting of switches simul- 
taneously operating with a period 7, and advance and 
delay elements as shown in Fig. 3. The validity of this 
model can be demonstrated by purely physical reason- 
ing. Fig. 4(a) shows a switch S sampling an input f(é) 
with the period 1/3 to yield the sequence fo, fi, fo, etc., 
at its output. In Fig. 4(b) the same function f(t) is 
sampled by a combination of three switches S,, S:, and 
Ss, in parallel, each of them operating with the period T 
but with a time delay 7/3 with respect to one another. 
Output sequences corresponding to every switch are 
shown in the diagram and it is obvious that their sum is 
identically equal to the output of the switch S. It is 
always possible to simulate the effect of a delayed opera- 
tion of a switch by advancing the input function, sam- 
pling it without a delay, and then delaying appropriately 
the output sequence in order to place it in the correct 
time in relation to the rest of the system. Therefore, the 
operation of the switch, say S:, can be simulated by 
placing advance, e*7/*, and delay, e—*7*, elements at the 
input and output of a switch which closes simultane- 
ously with the switch S. In this way all delayed switches 
can be replaced by switches without any time delay but 
with appropriate delay and advance elements, and 
therefore any switch operating with a period T'/n can be 
replaced by a system of switches operating simultane- 
ously with the period 7, as shown in Fig. 3. 

Assuming that a sampled-data system to be analyzed 
contains switches operating with periods 7), To, -- - , 
T,, it is always possible to find a number JT and a set of 
integers pi, po, - ++, pn such that 7) =7T/pi, T2=T/po, 

-++,7,=T/pn, provided that the ratios of 7; to T> to 
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Fig. 4—Signals in switch decomposition procedure. 


- T,, are expressed by rational numbers. This limita- 
tion is not serious since we can always approximate an 
irrational number by a rational one. It follows from the 
previous discussion that the system under analysis can 
be replaced by an equivalent mathematical model, with 
all switches operating with the same period 7, thus per- 
mitting the use of periodic properties of Z transforms in 
the formulation of the loop equations. The “switch de- 
composition technique” described above is the key to 
the problem under discussion. 


General Case 


The feedback system shown in Fig. 5(a) will be con- 
sidered. It contains two switches operating with periods 

=T/pand T;=T/q where p and g can be any integers 
(excluding of course zero). Since, as will appear later, 
the presence of additional switches would only increase 
the amount of algebra necessary without affecting the 
method of procedure, this case deals with all the basic 
points encountered in the general problem. The first 
step in the procedure is to apply the “switch decomposi- 
tion technique,” thus obtaining the mathematical 
model shown in Fig. 5(b). The following two equations 
can be formulated: 


E(s) = R(s) — B(s) (16) 


B(s) = aN + Z[G(s)e? 4] H(s)es7/0 
Z[G(s)e2*T/9] A (s)e28T 1a 4 
| De eh aaa 
Z[E(s)e??]{ Z[G(s)e*7/”] H(s) 
Z[G(s)e~8T/P eT 4] H(s)e rece eke 
+ Z[G(s)e*? Pe sT 14] AT(5)e-(g-Y eT a} plata k 
Z[E(s)e@-Y eT ln] { Z (G (s)e—@—D eT 1] (5) 
Z[G(s)e— ODT Imes Ta] A (seo 4 +. 
Al 


G(s)e~@-Y)sT Ipg(a-1)8T Ja] FI(s)e— 2/4} (17) 


T | 
(b) 


Fig. 5—Multirate feedback system; sampling rates related by 
any rational number. 


Applying Z, transformation to both sides of (16) and 
(17), the following relations are obtained: 


E(2p) = R(zp) — BZ») (18) 
B(Zp) = E(2)gi(zp) + Z[E(s)e7!?]g0(2p) + - 
ef ae |S ew tr eae) (19) 


where 
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q—t 


gilt») = De Z[G) Fe] Zp (see 9] 


glee) =O Z(G eet a] Z, [21(s) bt 


A qt 
= Dd. Z[G(s)e— @-YsTlngksT /a] 


k=0 


(Zp) Zi\H(s)eeF ia] (20) 


Eliminating B(z») from (18) and (19), the following re- 
lation is obtained 


E(2) — E(z)gi(zp) — Z[E(s)e*?!?]g0(zp) — 
m= Wintejce-bitisly (5) = = R(zp) (21) 
where E(z,), E(z), Z[E(s)e7/?], Z|[E(s)€@-Ds7/? | 


are unknown terms. To solve (21) for E(z,) use is made 
of the identity derived in Appendix I, namely: 


1 P 
Z[E(s)eeT/?] = 7 Spt >, €28a(2—n) IP F(z e128 (v—n)/P)_ (22) 


n=1 


The above identity permits the expansion of all the un- 
known terms in (22) into functions of z,. Performing the 
indicated expansion and rearranging the terms the fol- 
lowing relations are obtained: 


E(Zp)fi(p) + E(zpe??*!?)fo(zp) + - > 
t E(Zpe4?—)!?) f,(Zp) =s R(Zp) 
or more concisely 


DE efr(k- DIP) (Zp) = R(Zp) (23) 
ALS E: p 
where fila) = 1+ 2 ao 81(Zp)Zp' + 
toe: 
falep) = — DD gulep)[zpei er! ]> 
Pp Oat 
ee? 
frl%p) = ep re gi(Zp) OP cpccmilt Bhs (24) 
l=1 


Eq. (23) is a functional equation with the term 
E(Zpei2*%-Dip), k =1, 2, -, p representing p unknown 
quantities. Additional p—1 equations with the same un- 
knowns are readily formed by substituting z,€??"/, 
- , gpei2"@-VIP for gz» in (23). In this way, the 
following set of p simultaneous equations with p un- 
knowns is obtained: 


Zest, aeie 


Sy E(zye!?* 010) f(a) = R(Zp) 


k=1 


_: 2 
» E(zpe??™*!?) fy(Zpei?*/?) = R(zpe?"/?), 
k= 


; 53 E(z,e12" F229) I?) ( f, (gpet** P—))/?) = R(Zpef?* GV?) (25) 


; 
k=1 
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Eq. (25) is solved then, algebraically, for the unknown 
function E(z,) and the output sequence is given by the 
relation: 


C(Zp) = E(2Zp)G(Zp) . 


It is possible to replace loop equations (16) and (17) by 
equations where the transform of the output rather than 
the error signal, appears as the unknown. Applying the 
Z, transformation to these new equations and making 
use of the identity (22), C(z,) is obtained as the solution 
of the functional equation. Since, in general, this in- 
volves g simultaneous equations with g unknowns, this 
procedure is advantageous whenever g is less than p. It 
is observed from (20) that in order to evaluate auxiliary 
functions g,(z,) it is necessary to find expressions of the 
form Z| F(s)e*7] and Z,|F(s)e #7]. The problem of find- 
ing Z[F(s)e8?] in a closed form has already been 
treated in the literature’ in connection with plant 
delay lags. A brief resume of the procedure to find 
Z| F(s)e7 |, and its extension to the case of Z,|F(s) 7], 
is given in Appendix II. 


(26) 


ik 
Tog 


Fig. 6—Illustrative example. 


Illustrative Example 


The feedback system shown in Fig. 6 is to be consid- 
ered as an example illustrating some of the points dis- 
cussed in this section. In this case 7; =1/2 sec, T2=1/3 
sec, G(s) =1/(s+1), H(s) =1/s, and R(s) =1/s. To take 
care of the effect of the small but finite time 7 during 
which switches are assumed to be closed, G(s) and H(s) 
are defined here as the Laplace transforms of the re- 
sponses to impulses of strength 7 rather than unity. 
Since T; and T. can be written as T/2 and T7/3, re- 
spectively, with T=1 sec, it follows that p=2 and 
q=3. Eq. (25) reduces, therefore, to 


E(@2)fi(g2) + E(—22)fo(z2) = R(z2) 


(27) 
E(22)f2(—22) + E(—22)fi(—22) = R(—22) 
and E(ze) is 
R(z2)fi(—22) — R(—22)f2(22) 
fo file)fi(—22) — fo(—22)fo(z2) se 
From (24) 
filze) = 1 + 3g1(22) + 3228222) 
and 
fol%s) = 381(22) — 3%2g2(22) (29) 
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where 
gi(z2) = Z[G(s)]Z2[H(s)] + Z[G(s)e?/?]Z2[H (se?! | 
+ Z[G(s)e27!3]Zo[ (5) e287 /3] 
go(z2) = Z[G(s)e*?/?|Z2[H(s) | 
+ Z[G(s)e*7/8|Zo[H(s)e87 3 
+ Z[G(s)e7!*]Z.[H(s)e*? 8]. (30) 


Each term in the above relations can be evaluated mak- 
ing use of Appendix II. Replacing all z by 2? finally ob- 
tains the following expressions for gi(z2) and go(z2): 
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sampling period 7. Without detracting from the gen- 
erality of the case under discussion, the feedback system 
shown in Fig. 7(a) is to be considered. This system con- 
tains only three switches, operating with the sampling 
periods T, T/p, and T/q, p and g being any integers not 
equal to zero. Applying the “switch decomposition” 
technique, the equivalent system shown in Fig. 7(b) is 
obtained. Since the error signal e(t) is the input to the 
switch which operates with the sampling period T, the 
loop equations are written in terms of F(z), as the un- 
known, and they read 


ak Mil a R(@) = RG) - BE) (33) 
gi(22) = 
(1 = romrarevme yn ell == Zy') and 
and B(z) = E(z){Gi(2)[G2(2)H@) + --- 
(as (e-7/6 4 e-T/2)g.—-2 1. ¢-5T/6g,—3 31) es Z[Go(s)e@- 87/9] Z[ H (5) e~ 2-8/4] ] Eure: 
§2\22) = } 
(Laem cia yin es-*) + Z[Gi(s)e@-YsT/2][ Z[Go(s)e- DTP] T(z) + ++ 
therefore, noting that R(z:.) =1/(1—227}), + Z[Go(s)e— @-D sT/P¢(2-D) 87/9] Z] A (s)e— (a-YsT/q] | } ‘ (34) 
il + i ea Z il fe 52/16) Let zg Pee, Ts 1) 
E(2s) = siete 8 ENE eae oe (32) 
24+ (-1 — eH eT! + ePT/3 + e-5T/6 4. ¢-T/6) 75-2 4 e-lz.—4 
To obtain C(z2) use can be made of (27) which reads __ or more concisely 
in this case, C(z2) = G(z2) E(ze). elk 
B(z) = E(2) D1 Dd) Z[Gi(s)e D477 ]Z 
SystEMS WHERE AT LEAST ONE SWITCH m=1 k=1 
OPERATES WITH PERIOD T - [Go(s)e— “D8 lpg RV) 87 1a] Z| H(s)e— Vs TI], (35) 


In this section a special case of a feedback system, 
containing at least one switch operating with the sam- 
pling period T, is considered. It is assumed that other 


E(z) = 


Eliminating B(z) from (33) and (35) the following ex- 
pression for E(z) is obtained 


R(z) 


(36) 


p q 
1+ Da Ny Z[Gi(s)e YT /P]Z[Go(s) es TPE HV) 87/4] Z| FT (s)e—F-VeT/a] 


m=1 k=1 


switches in the system operate with periods 7/f,, 
T/po,--+, T/pn, where p1, po,+ ++, Pn can be any 
integers not equal to zero. The property of the case dis- 
cussed here is that it is possible to write the loop equa- 
tions so that the solution of the problem does not lead 


Finally, since 
C(s) = G2(s)Zp[Gi(s) E(2)] = Ga(s)Gi(zp) E(Zp?), (37) 


and therefore 


to functional equations. The loop equations ought to be C(Zp) = Go(2p)Gi(Zp) E(2>”) (38) 
set up in terms of the unknown variable representing ve 
the signal at the input to the switch operating with the it follows that 
R(Zp?) Gi (Zp) Go(z 
C(Zp) ba ( ‘p ) 1( ) a( p) 


Pp q 
1-2 es Z[Gi(s)e™ 57 !?]Z[Go(s)e— MDs Plag&—1) 57/4] Z| A (5) e- VT Ia] 


k=1 k=1 


(39) 
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Fig. 7—Multirate feedback system; at least one switch operating with period T. 


where z,” is to replace z in the terms of the denominator. 


CONCLUSION 


The analytical technique described here makes it 
possible to study the behavior of multirate sampled- 
data systems. It has been shown that, in principle, 
it is always possible to obtain the closed-form Z trans- 
form of the output of such systems and, therefore, 
standard tools of sampled-data theory could be applied 
to deal with the problems of stability and performance. 
The amount of computation necessary to obtain the 
closed-form Z transforms of the output depends in the 
feedback cases on the actual values of the integers » and 
g (assuming that switches with periods T/p and T//¢ are 
present). It is obvious that the higher their values, the 
more computational labor necessary. This becomes es- 
pecially serious in the general case, where it is shown 
that in order to obtain the Z transform of the output 
(p or g, whichever is smaller), simultaneous algebraic 
equations have to be solved. Computational labor could 
be reduced by suitable computers or the use of approxi- 
mations. Thus, it seems reasonable that it is possible to 
obtain some estimate of the system performance by ap- 
proximating sampling periods, which lead to high values 
of integers p and g, and by sampling periods for which 
the values of p and g are small. In cases when p is large 
and g is small or vice versa, it seems possible to ap- 
proximate the ratio p/q (or q/p) by an integer n, thus 
reducing the system to a case which is easily computa- 
ble. In this case (sampling periods related by an integer) 
the amount of computation necessary is independent of 

the actual values of the sampling periods. 
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To prove that 


p 
Spt se ela(p—n) 2x /p Fi] 2 i (P—n) 20 /P | 


1 
Z[E(s)ee?/2?] = — 
p n=1 


where g and # are any positive integer, p not equal to 
zero. 

Proof 

Let e(t) = £—1E(s). Then 


ve) 


A 
Z[E(s)et?2!7] = DY e(RT + g/pT)z+ 


k=0 


: D 2(s+in— (stin2e/T)Ta/ 
a aS, 5 v— stjn2r alp 
Tee. ORS 


by Poisson summation rule. Since e*?/? =g,%, then 


ico} 


E $fa/p) = = pi pet jn2rq/ 
Z\E(s)e*Fa!?) =. — z r 2 Bs +52 am) ¢ 4 


Since, in general, 


eke pa ee Os 


n=—oo n=kp n=kp+1 n=kp+p—1 


where R=—0o,--+—1, 0, 1,-+- 0, and p equals 
any integer not equal to zero, therefore, 
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or) 9) , 
Sy (s + jk =) eikade 
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L\E SLD) | Bs ee 
[eG )erte*] me — sgt De 


co i : 
a B| s+ +i(® - =| ee <| ei (ktl/p)gart ... 
a2) Bs a i(e ic i) 2r f elle (o— vine ; 


Since, e*22"=1, the above relation can be rewritten as 


E a(+0% 


k=—0o 


it 
Z| E(s)e?a/? | = — ate 
[E(s)e*7*/?] of 


= 1 
+ etarm £ ye HE +i(k+—) 2" Fl+- 


k=—oo 


= 
+ comune © 2) zs +o (i ob 4 QT = It ; 


or more concisely, 


nas at ef (a/p) (p—n) 2x 
p- 
eS). 
= ko ee | 
os stj(k+ a) a 


Z[ Bet] = 


fh 


To complete the proof it is sufficient to observe that 
since 


Be) * See D)ar= 2 ¥ 2s +m 2), 


k=—oo 
therefore, 
p< Pia ") | 
pee ieee Y Sect 
F als ti(at p ? hs 


= E[zpeie-n/?) 2], 


APPENDIX II 


In Table I, reproduced from Bergen and Ragazzini,‘ 
to find Z[F(s)e*?] apply the following rules: 


Dire b> set ya = -A,; 
therefore, 


Z|F(s)e8?| = Z[F(s)e47*]. 


2) a> 1; set a=m+A 


where m is a positive integer smaller than a. 


Therefore, 
Z|F(s)ew*?] = Z[F(s)em*TeAT?] = 2Z[F(s)e7*] . 


3) =n Ss a< 0 
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TABLE I 
Z-TRANSFORM PAIRS FOR FUNCTIONS OF THE Form F(s)e478,0<A <1 


F(s)eA7s Z [F(s)eA7s ] 
i ATs 1 
§ 1—z-1 
hs AT+Tz71(1—A) 
—- ATs 
S$ (1—z—1)? 
1 @ @AT 
——— ean ps es Re 
sta 1—eaTz1 
a rem itr: aAT 1» 1(€ aAT _ ¢ aT) 
s(s+a) (1—23)(1—e*7?z) 
where 7 is a positive integer; set a= —h-+A. 
Therefore, 


Z[B(s)e*™ | = Z[F(s)e meth?) =e *Z| Ps) At: 


In all cases the problem of finding Z[F(s)e7] reduces 
to finding the expression for Z| F(s)e47*] from Table I. 

To find Z,|F(s)e%*7 |-frst note that it can be rewritten 
as Z,|F(s)e*?s7/?] and then apply the following rules: 


4A) 0 Sap =< 15" set ap= A, 


and therefore 


Z,[F(s)erPetl?] = Zp|F(s)e4#7/?]. 


3) ape; set apf=m-+-A 


where m is a positive integer smaller than ap. 
Therefore, 


Z[F(s)er?2T!? | ES Zp[F (sem? Phe? lp] ae t2"Zy|F(s)et*t lz], 
6), =n < ap <0, 
where 1 is a positive integer. 
Let ab = —n-+A and therefore 
Zp[F(s)err*? 9] = Zp[F(s)em? veto 2] 
2, *Z pl F(s) Att /0]., 
Thus the problem of finding Z,[F(s)e*7] reduces to 
finding the expression for Z,[F(s)e4:7/7]. To do this it is 
sufficient to obtain from Table I, the expression for 


Z|[F(s)e47*| and then to replace every T in it by T/p, 
and every 2 by 2p. 
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Statistical Design and Analysis of Multiply- 


Instrumented Control Systems” 
ROBERT M. STEWARTT 


Summary—It is the purpose of this paper to show how Wiener’s 
linear least-square filter theory may be applied to some common 
types of multiply-instrumented control systems. By multiply-instru- 
mented control systems are meant those in which more than one 
sensing instrument is used. 


GENERAL SOLUTION AND EXAMPLE OF Two- 
Loop SERVODESIGN PROBLEM 


| Panes 1 depicts a fairly general form of two-loop 
servo in which D isa disturbance; €, and e€: are errors 

in two controlled variables; N; and Np» are noises 
Or errors in measuring & and €, respectively; Y;, and 
Y;, are two transfer functions characteristic of the 
“load” being controlled and are ordinarily specified; 
while Yc, and Y¢, are the controller transfer functions 
which are to be designed. If D, Ni, and N2 are mutually 
independent stationary random variables, then by 
straightforward servo and harmonic-analysis tech- 
niques it may be shown that the spectral density of €; is 


Bie Aa sii | Vz,|*h0} + 


eee 


Gail 


X,|?oy, 


where 
Vy Ve, palte(1) 
See f hal Cat okt ee 


Doge 


and 
a Vel a, 
i cP Bc, 


The mean-square error €? is given by the integral 
over all frequencies of ®., and the forms of Xi and X» 
(and, hence, Yo, and Ye,) which minimize. It may, in 
principle, be found for any ®p, @y,, and ®y, by using 
the technique outlined by Wiener! since (1) corre- 
sponds exactly to the two-input filtering problem. An 
example which has been worked out is one in which 
acceleration and position measurements are used simul- 
taneously to minimize deviation in position of a load of 
primarily inertia character.” Then, in transform notation, 


* Manuscript was received by the PGAC, November 10, 1956. 
The work was supported by the Army Ordnance Corps, ORDCIT 
Project, under Contract No. DA-04-495-ORD 18. Presented at 
WESCON, Los Angeles, Calif., August, 1956. 

t+ Jet Propulsion Lab., Calif. Inst. Tech., Pasadena, Calif. 
1N. Wiener, “Extrapolation, Interpolation, and Smoothing of 

Stationary Time Series with Engineering Applications,” John Wiley 
& Sons, New York, N. Y.; 1949. _ i 

2 This is a possible representation of the problem of continuously 
guiding a missile by means of a combination of inertial (accelerome- 
_ter) and radar data. 


Fig. 1—Dual-loop control system. 


1 

p 
For the case in which both instrument errors have flat 
spectral distributions over the range of frequencies of 
interest and 


Yan 


¢ 
aintan=. me 


where @p is a constant, the optimum controller transfer 
functions turn out to be 


iol Dy, 
Le -)- | 
2 
fies a ary @ 
ly ge- 8-2] 
where eke) 


T°>/ bp Py, Py, 


®y, by, 
| ar4/ — 14 =U) B 
op ®y, ( 2) 


2+ 7? — 
Py, 


2by,Py, 


“if 
op —) 
pe 2 ee 
(2 ala a ( a by, 
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It should be clear from this simple example, however, 
that the utility of the method is limited somewhat by 
its complexity to the simplest problems. Thus, current 
and future research effort, applied to methods of sim- 
plifying this design technique even at the expense of 
slightly less precise results, will be extremely useful. 


DISTORTIONLESS MEASUREMENTS 


A technique which has been used frequently in sys- 
tems, such as the one discussed in the first section, is that 
of distortionless or “complimentary” measurements. 
(See for example Ryerson* where, however, only the 
optimum “unrealizable” solution is given.) In this case, 
there is no contribution to error associated with the 
signal being measured, but only from the instrument 
errors. 

A typical example is shown in Fig. 2 where S* is a 
continuous estimate of the signal S. 


S+N, 


Fig. 2—Typical distortionless system. 


This mixing circuit may be shown to be the optimum 
(in the Wiener sense) distortionless system when the 
spectral distributions ®y, and ®y, of Ni and WN; are both 
flat. Then, the optimum time constant is given by: 


Py, ye) 
Tmo = 
a laa 


and the corresponding mean-square error is: }. (3) 


(¢2)min = —— [y,4y,]*”? 

A legitimate comparison with a single-input system 
often is complicated by various special considerations; 
but if it is assumed that the spectral distribution of 
dS/dt, (®g), is flat, then the ratio of minimum mean- 
square error with input one to that with two compli- 
mentary inputs for the case above is 


oi Gg 1/2 
“(27 
on) Py, 
The ratio on the right is essentially the signal-to-noise 
ratio of the derivative channel. 

Fig. 3 shows two equivalent general methods for 
mechanizing distortionless dual-input systems with a 
single frequency-sensitive element. 

The error-spectral density in any case is given by 


3 J. L. Ryerson, “Trajectory precision requirements for automatic 
landing,” IRE Trans., vol. ANE-2, pp. 4-10; March, 1955. 
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Fig. 3—General distortionless dual-input systems. 


= | Vi[’éy, + |1— Yi |?4y, (5) 


and the frequency response functions at the output 
relative to each input are related by 


Vi(jw) + Yo(jo) = (6) 


Thus Wiener’s solution for the single-input case may 
be applied to (5) to determine the best form of Y, and 
then, if desired, the corresponding Y2 can be determined 
from (6). 


An APPLICATION OF WIENER’S MULTIPLE- 
INPUT FILTERING THEORY TO 
DIVERSITY SYSTEMS 


Another technique of considerable interest in con- 
temporary communication, control, and instrumenta- 
tion systems is illustrated in Fig. 4. 

Wiener’s general analysis of such problems shows 
that, 1045, 00VG; Van , VN, are mutually independent 
stationary random time series, the impulsive responses 
yi(r) of the linear filters which minimize the mean- 
square error in estimating the “signal” S by S* satisfy 
the following simultaneous set of integral equations: 
G=1, 2, a - 2) 


os(t) = [ Xlec —s)] 2 9 


+ CG — c)loo)t do forte > Ona 


where the ¢ terms are autocorrelation functions of the 
signal and various noises. 
These equations are equivalent to: 


boo n 
alr) — fi {lost = 0] 90) 
sep Reteie e)|yi(o)} do = h(n) (8) 


with the understanding that 


yi(c) = 0, tObar= 0; 


(9) 
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Fig. 4—Diversity system. 


(the “realizability” condition) and 


hj(ry = 0, ies 5 0: (10) 


Fourier-transforming both sides of (8) then gives 


#3(0)[1— Ve) | = ey) Mle) = He) AY 
t=1 

where the ®’s are spectral densities of signal and of the 
various noises. The Y’s are frequency-response func- 
tions of the optimum filters having [from (9)] no poles 
in the lower half-plane of w, and the H’s having [from 
(10) | no poles in the upper half-plane.4 

A frequently encountered situation is that in which 
the noises in the various channels are statistically simi- 
lar, except possibly for amplitude. Then, 


Py, = K;®y,. (12) 
Substituting (12) into (11) then gives: 
#3| 1 —-> r.| — K7y,V; = Hy. (13) 


Subtracting each side of this equation for arbitrary j 
from that for 7 =1, gives 


by,[K;?V; — KV,] = Hi — Wj. (14) 
Using Wiener’s spectral-factorization theorem, 
Py,(w) = we) = yy (15) 


where ¥* and y~ have poles and zeroes only in the upper 
and lower half-planes, respectively. Thus, 


Hie Ai 
lip 
The left side of (16) can then have poles only in the 


upper half-plane, while the right side has poles only in 
the lower half-plane; hence, 


¥+[K?Y; — KY,] = constant. 


¥+[K7?V; — KYY,| = (16) 


(17) 


4 Corresponding poles of s=jw, rather than w, would be in the 
right and left half-planes, respectively. 
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As wo, all | Y;| —0 since they have been assumed to 
correspond to well-behaved weighting functions y,(c) 
and must, therefore, be Fourier-transformable. Thus, 


v+|K?V; — KLYi] 0 (18) 
and, in general, 
Ki Vn Kav (19) 
K? 
Vile) = (55) 7. (20) 
IK 


Eq. (20) indicates that the optimum system may be 
mechanized with just one frequency-sensitive element 
as shown in Fig. 5. This implies the somewhat surprising 
result that the optimum bandwidth for each channel is 
identical, even though the signal-to-noise ratios may be 
quite different. 


S+N, 3 
S+Np Kae 
>) 
K>° 
* 
Ss 
° Y, () |\—» 
S+N, K,2 : 
Kye 


Fig. 5—Possible mechanization. 


The best form of Y; may be determined as follows: 
substitution of (20) into (13) gives 


ted 
bg E — KY, >> =| —K/y,Vi=H, (21) 
Or, 
by — [bs + Sy,/|V = Ai (22) 
where 
® 
Oat aman (23) 
ae 
and 
Vi =| K? 53 : Y (24) 
| Nee an 1 Kz 1 


Eq. (22) may be solved for Yi’ by the method applic- 
able to a single-input filter, and then Y; may be found 
from (24). 
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A Time Domain Synthesis for Optimum Extrapolators" 
CARL W. STEEG, JR.+ 


Summary—A direct method is presented for the solution of the 
integral equation for the optimum predictor in terms of the solution to 
the integral equation for the optimum filter. This method is a means 
for circumventing the practical difficulties encountered in the actual 
design of optimum predictors based upon the techniques derived by 
Wiener. The synthesis procedure is applied in order to simplify the 
design of extrapolators for use when the prediction interval is a non- 
negative function of time in contrast to the more conventional situa- 
tion where prediction is made for a continuously varying instant that 
is always a fixed number of seconds in the future. The specific exam- 
ple utilized is prediction for a fixed instant that is a continuously de- 
creasing number of seconds in the future. The discussion includes an 
explanation of a procedure for avoiding the solution of integral equa- 
tions in the synthesis of optimum extrapolators. 


I. INTRODUCTION 


NE important recent change in communications 
() theory and control-systems engineering is the ex- 

tension of interest from primarily static perform- 
ance characteristics to those of a system excited by time 
functions including actual typical random inputs. Tech- 
niques for the design and synthesis of control and com- 
munications systems also have been extended to include 
consideration of both time-domain and frequency- 
domain characteristics. This paper presents a typical 
problem in automatic control theory that represents 
the trend away from frequency-domain methods. 

The stationary statistical optimization techniques de- 
veloped by Wiener! and others?:* and subsequently ex- 
tended by Booton‘ to time-varying systems not only 
have made valuable contributions to communication 
theory but also have been a stimulus to the application 
of statistical concepts in the design of other types of 
systems. In these fields, design often requires that the 
future value of a signal be predicted. Either of two dis- 
tinct situations frequently may arise: 1) prediction for 
a continuously varying instant that is always a fixed 
number of seconds in the future, 2) prediction for a 
specific fixed instant that is a continuously decreasing 


* Manuscript received by the PGAC, December 12, 1956. The 
work was supported in part by the U. S. Air Force under contract 
no. 33(616)-3398, monitored by the Aeronautical Res. Lab., Wright 
Dey. Ctr., and in part by the Bureau of Ordnance, Dept. Navy, 
under Contract No. NOrd 16594, both administered through the 
Div. of Sponsored Res., Mass. Inst. Tech. This paper was presented 
at WESCON, Los Angeles, Calif., August, 1956. 

{ Radio Corp. of America, Waltham, Mass. 

1N. Wiener, “Extrapolation, Interpolation, and Smoothing of 
Stationary Time Series With Engineering Applications,” The Tech- 
nology Press, Cambridge, Mass.; 1949. 

2S. O. Rice, “Mathematical analysis of random noise,” Bell Sys. 
Tech, J., vol. 23, pp. 282-332; July, 1944 and vol. 24, pp. 46-156; 
January, 1945. 

3H. M. James, N. B. Nichols, and R. S. Phillips, “Theory of 
Servomechanisms,” M.I.T. Rad. Lab. Ser., McGraw-Hill Book Co., 
Inc., New York, N. Y., vol. 25, 1947. 

*R. C. Booton, Jr., “An optimization theory for time-varying 
linear systems with nonstationary statistical inputs,” Proc. IRE, 
vol. 40, pp. 977-981; August, 1952. 


number of seconds in the future. The first situation, 
which is the classical stationary prediction problem, can 
be treated by the techniques derived by Wiener.’ The 
second situation demands nonstationary theory and is 
suitable for solution by the theory given by Booton.* In 
each instance, an essential step requires the solution of 
an integral equation that involves the impulse response 
of the optimum system. The usual procedure determines 
the Fourier transform of the impulse response of the 
physically realizable system through the application of 
frequency-domain techniques. Although Bode and 
Shannon’ have given a thorough discussion of the sta- 
tionary prediction problem, a similar description of the 
nonstationary prediction problem has been lacking. 
These two prediction problems are merely special cases 
of the general nonstationary problem of predicting over 
a variable interval. 

The application of statistical optimization theory to a 
system is suitable whenever the input to the system 
consists of a random signal plus a random noise and if 
the system is to be designed to produce a response that 
approximates a desired response. The signal and noise 
are assumed to be nonstationary in the sense that each 
are derived by means of a linear time-varying integro- 
differential operation on white noise. Although this re- 
quirement appears to restrict the applicability of the 
results, it serves both to simplify the mathematical 
equations and to facilitate the physical description. 
Furthermore, the only physical systems that are amen- 
able to treatment by existing optimization theory belong 
to this category of nonstationary systems. The response 
of the actual system is compared with the result of a 
desired, time-varying, linear operation upon the signal 
component of the input, and the resulting difference is 
called the error. Because time-varying operations are 
admissible, the entire system must be considered in 
terms of nonstationary theory,‘ in spite of the assump- 
tion that both signal and noise originate as white noise. 
The optimization criterion is taken as the ensemble- 
averaged mean-square error, and the characteristics of 
the physically realizable system are chosen in such a 
way that this mean-square error is minimized. An as- 
sumption inherent in the use of mean-square error is 
that the system is linear. A linear system belongs to one 
of two broad categories: 1) systems described by finite- 
order differential operations and 2) systems not so de-_ 
scribed. The first class of systems is predominant in 
automatic control theory. The second class is seldom en-_ 
countered except in abstract theory of communica-_ 


_ °H. W. Bode and C. E. Shannon, “A simplified derivation of : 
linear least square smoothing and prediction theory,” Proc. IRE, 
vol. 35, pp. 417-425; April, 1950. | 


1957 


tions. Because the input to the system contains noise, 
a filtering problem exists even though the desired opera- 
tion may be merely reproduction of the input. If the 
result of the desired operation depends upon future 
values of the input, the fundamental filtering operation 
must be supplemented by a prediction or, more pre- 
cisely, an extrapolation. The system described is shown 
in Fig. 1, where the fundamental dependent variables 
are defined. 


v(t) 
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Fig. 1—Block diagram of optimum extrapolator. 


Prediction problems usually are classified according 
to whether nonstationary or stationary statistical opti- 
mization theory is needed for the solution. This paper 
proves the equivalence of the two prediction problems 
and derives a procedure for determining the optimum 
predictors. The procedure requires knowledge only of 
the optimum filter and its impulse response. For this 
reason, the system designed is called an optimum extrap- 
olator to distinguish it from the conventional optimum 
predictor. The extrapolator itself is a device that oper- 
ates upon signals present in the optimum filter to obtain 
the desired predicted value of the signal component of 
the input. First, the techniques are established on a 
rigorous mathematical basis in the time domain. Next, 
an interpretation of the significance of each signal pres- 
ent in the predictor is given, together with a typical 
example that includes the results of an analog-com- 
puter study. In conclusion, the potential applications 
and usefulness of the theory are discussed. 


II. THe MATHEMATICAL THEORY OF 
OPTIMUM PREDICTORS 
If x4(¢) denotes a statistical function, the autocorrela- 
tion function yu is defined by the ensemble average 
yaa(l, te) = xa (t1) xa (te) (1) 


where the bar denotes the averaging operation. If the 
signal is stationary, the autocorrelation function de- 
pends only on the difference between ¢, and ¢, and is 
written in the more familiar form 


gaa(T) = wa(h)xa(tr +7). ; (2) 


Similarly, the cross-correlation function yas between 
two statistical functions x4 and xz is defined as 


2 yanh, bo), = wa(ts) a (te). : (3) 
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Booton! has established that minimization of the 
mean-square error requires that u(t, fe), the impulse 
response of the system defined by the superposition in- 
tegral, 


xtr(t) = fu to) x1(to) dte (4) 


u(t, te) = 0, iy 2 Hl 


be determined by the solution to the integral 


t 
yrp(ti, t) = f u(t, te)yrr(t1, t2)dte, See EM 


et) 


where yzp is the cross correlation of the input x; and the 
desired response *p and yz is the autocorrelation func- 
tion of the input. To provide consistency for the current 
development, some slight modifications have been made 
in the notation used by Booton.! 

The integral (5) can be reduced to simpler terms ac- 
cording to the characteristics of the input and of the 
desired response. For example, the function of an opti- 
mum filter is to reproduce the signal component of the 
input; that is, the desired response is the signal. Thus, 
if the signal and noise components of the input are un- 
correlated, the impulse response of the optimum filter 
g(t, t2) is determined by the solution to the integral (5) 
with yrp equal to yss, 


t3 
a g(ts, te)yrr(t, t2)dte, ti < te | 
ie (6) 


t3 < lo. 


vss(t1, ts) 


g(ts, te) = 0 


Similarly, the desired response of an optimum predictor 
is equal to the signal evaluated at some future time 
t+a, and the impulse response of the optimum predictor 
h(t, t2) is obtained by the solution to the integral 


t 
vss(t1, t+ a) = ir h(t, te) yir(ti, ty) dte, tt, (7) 


h(t, t2) = 0, fs, 


where a@ can be any nonnegative function of ¢. If the sub- 
stitution t3=t-++a is made, (7) becomes 


t3—a@ 
ae tp Nilo TOMA ARON EE: 
ai ty <a t3 es (8) 


h(ts ae (35) to) =0(, tz < le =e Qa. 


Because the latter part of (8) holds, the upper limit on 
the integral can be written ¢; rather than ¢#;—a@. Thus, 


43 
yss(ti, ts) = i} h(ts — a, te)yrr(t1, te)dte, ti <<ts— a. (9) 


Design procedures! for determining the optimum 
filter provide the necessary means for ascertaining the 
function g(ts, tz). The following development describes 
a method for the determination of the optimum pre- 
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dictor when only the impulse response of the optimum 
filter is known. The response xz of the optimum filter is 
given by 


xr(t) => ic t1)¢1(t1)dty (10) 


where g(t, 1) satisfies (6). Similarly, the response of the 
optimum predictor can be written 


Xxep(t) = fn ty)ar(t1)dty (11) 


where h(t, t:) satisfies (7). 

Although this equation often is used as a basis for the 
design of optimum predictors, a slightly different con- 
figuration for the system offers advantages in the inter- 
pretation of the duality between the stationary and 
nonstationary prediction problems. The configuration 
preferred for this purpose is shown in Fig. 1 where the 
straight transmission corresponding to g(t, f1) is replaced 
by an equivalent unity-feedback transmission involving 
y(t, t:), the impulse response of the forward path. The 
configuration represents a significant departure from 
the techniques used in the references. Reasons for the 
choice of this configuration are emphasized in the de- 
velopment of Appendix I. The following definitions, 
indicated in Fig. 1, are convenient: 


xn(t) = x1(t) =— xp(t) (12) 
and 
xR(t) => ic ty) Xp (ti) dty 
= fx Ale) = eats lan. 4) 


To be physically realizable, the impulse response of the 
forward path y(t, ¢,) is zero for t, >t. 

The remainder of this section is concerned with es- 
tablishing the validity that, under the condition that 
the noise component of the input is white, 


<rp(t) = i; y(t + a, t)en(tr)dty (14) 


—7o 


relates the response of the optimum predictor to the 

error signal of the optimum filter and to the impulse re- 

sponse of the forward path of the optimum filter. 
Because 


tta 
vat ta) = [ y(t + a, t)en(tr)dts 


menic°} 


t 
= { y(t + a, t)an(ti)dty 


= OO) 


t+a 
te i y(t + a, t)«n(tr)dtr, (15) 
t 


the expression (14), for the response of the optimum 
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predictor, indicates that the predictor response approxi- 
mates the future value of the filter response by the 
omission of the integral 


t+a 
f y(t + a, t1)xn(t1)dt. (16) 


This integral requires values of the error in the future; 
that is, for 4; in the interval between ¢ and t-+a. Because 
such values are not known, the integral must be neg- 
lected in the determination of future values of the in- 
put. This conclusion can be reached on a more rigorous 
basis by the argument given in Appendix I. 


II]. THe Trme-DoMAIN SYNTHESIS PROCEDURE 
The function of the optimum predictor is derived in 
the preceding section, where the output of the predictor 
is proved to be given by 


Xrp(t) = f y(t + a, ty) “n(t1)dty (17) 


68) 


where y(t, f;) is the impulse response of the forward path 
of an optimum filter and xz(¢) is the error signal in the 
filter. In an important class of optimum filters, the im- 
pulse response y(t, 4) can be written in the form 


n 


y(t, ty) = D2 es(t)b;(); 


j=0 


(18) 


that is, the impulse response, or kernel, is separable into 
the sum of the products of functions of f and 4 only. For 
this special, but dominant class, the optimum filter can 
be expressed in terms of a differential operator involving 
xr(t), xu(t) and their derivatives. The use of mean- 
square error as a criterion restricts the filter to a linear 
form and virtually all such forms are described by im- 
pulse responses that are separable. When the impulse 
response is separable, an expression equivalent to (17) 
can be established in order to simplify the design of the 
optimum predictor. The expression depends upon the 
impulse response of the optimum filter and involves 
only current values of xr(t), xz(t) and their derivatives 
in order to determine future values of the filter response. 
For this reason, the operation is called extrapolation, 
and the resultant system is termed an optimum extrap- 
olator. 

The equivalent expression for the integral (17) is ob- 
tained most easily from certain identities derived in the 
theory of ordinary differential equations.* Because these 
results are apt to be unfamiliar to engineers, the follow- 
ing discussion includes somewhat greater detail than 
might be required in establishing the desired equiva- 
lence. However, some facility in the manipulation of the | 
following equations is useful in applications such as the 
example given at the end of Section IV. 


°E. A. Coddington and N. Levinson, “Theory of Ordinary Dif- 
a Equations,” McGraw-Hill Book Co., Inc., New York, N. Y.; 


oy. 


If the coefficients of 


L(x) = So a) — = 0 (19) 
i=0 dt* 
are all finite, single-valued, and continuous throughout 
a prescribed time interval, the fundamental existence 
theorems for the solutions to ordinary differential equa- 
tions prove that there exists a unique solution x(t) such 
that x(¢) and its first »—1 derivatives assume a set of 
arbitrarily assigned values xo, x0, - + +, «9%, when 
t=to, and such that x(¢) may be developed as a Taylor’s 
series convergent in a certain interval around f.6-8 
Any linearly-independent set of 1 solutions to (19) is 
said to form a fundamental set or fundamental system 
of solutions. Although there is an infinite number of 
such sets, one particular fundamental set is of special 
importance because of its simplicity. Members of this 
particular fundamental set are designated as f;(t, to), 
7=0, 1,---, m—1, and these particular solutions to 
(19) are defined uniquely by the conditions 


ez bi, 


t=to 


d? 

= h(E; to) J = 0, i areal er uv (20) 
dt? 

where 6;; is the Kronecker delta, which assumes the 
value 1 when 7=j and is zero otherwise. Consequently, 


a unique solution to (19) is expressible as 


x(t) = x(to)fo(t, to) + x (to)filt, to) + °°: 
+ 4” (to) fn—1(t, to) 
a1 
= Do x (to)fi(t, to). (21) 
2=0 
Therefore, the derivation of the fundamental solu- 


tions f;(t, to) is sufficient for obtaining either the general 
solution to (19) or a particular solution satisfying pre- 
assigned conditions. The solution defined by (21) de- 
scribes the behavior that is usually called the transient. 

The condition expressed by (20) implies that the 
leading term of f;(t, to) is (¢—to)*/7! and that no terms in 
(t—to)**1, - - - , (¢—to)""! are present in the expansion. 
The usual method for the solution of (19) is to assume 
that 


ilfs(t, bo) = @ = 11+ 3 a(t = 1)", (22) 


j=n—t 


7 This statement of the existence theorems is valid only if é is not 
a singular point of the differential equation. Under the conditions 
outlined the only singular points that can occur within the prescribed 
interval are the zeros of the leading coefficient a,(¢). All other points 
are ordinary points. A singular point is defined as a pair of values of 
x and ¢ for which the solution to (29) either is discontinuous, is not 
unique, or does not exist. Because such conditions arise infrequently 
in control-systems applications, the assumption is made here that 
an(t) is different from zero throughout the entire time interval during 
which extrapolation operates. ; ; ‘ 
; 8 E. L. Ince, “Ordinary Differential Equations,” Dover Publica- 
tions, New York, N. Y.; 1944. z 
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to substitute the series in the differential equation, to 
arrange the resulting expression in ascending powers of 
(¢—to), and to equate the coefficients of like powers of 
(t—to) to zero. By this procedure, a set of simultaneous 
linear algebraic equations in the coefficients a;; is ob- 
tained. The relations between coefficients are called re- 
currence equations. The complete procedure can be cir- 
cumvented in many practical problems, such as con- 
stant-coefficient differential equations. 

The fundamental existence theorems for ordinary dif- 
ferential equations provide a basis for the prediction of 
many analytical functions. If T is a value of ¢ within an 
interval that contains no singular points, then, by (21) 


n—1 
x(t) = )) «(T)f,(t, T). (23) 
i=0 
Because relations between ¢ and T arébilateral, 
(ie " 
wT) = 2) eOWs(T, 0); (24) 
a—O 


that is, the value of x at any particular time is express- 
ible continuously as a time-weighted sum of values of x 
and its derivatives at any other time. Thus, in particu- 
lar, if T=t-+a, 


a(t + a) = si xD) (E + a, 2). (25) 


The definition of the fundamental set is such that, for a 
equal to zero, (25) reduces to the identity that x(¢) is 
equal to x(¢). As a increases from zero, increased weight 
is attached to the derivatives of x. 

For constant-coefficient equations which are invari- 
ant under translations in time, the members of the fun- 
damental set are functions of the difference between the 
two variables only. Hence, 


Filt, to) = falt — to) (26) 


for constant-coefficient systems, and the prediction 
formula becomes 


x(t + a) = yi x) (1) F(a). (27) 


If prediction is made over a fixed interval; that is, if the 
value of the function is desired for a continuously vary- 
ing instant that is always a fixed number of seconds in 
the future, the prediction formula is interpreted as ap- 
plying constant weights to current values of the deriva- 
tives of the function to be predicted. Consequently, the 
origin of the weighting functions encountered in the 
Wiener theory of stationary prediction is explained. 

Eq. (24) leads to interesting and useful identity 
among members of the fundamental set. Because the 
left-hand member does not depend on the current value 
of ¢, differentiation with respect to ¢ yields 
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7=0 g=0 
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= x(t) “jaT: t) + 2d, 2) 


+ fr—(F, t)ax™(Z). 


Because x(t) is derived from (19), the mth-order dif- 
ferential (19) and (28) must differ at most by a multi- 
plicative factor, designated as K(T, t). For example, 
comparison of the coefficients of the mth derivative re- 
sults in 


ad 
| Ha, ) + S07, | 


(28) 


ae t) ~~ ant) K(T, t). (29) 


The function K(t, t’) is called the unilateral Green’s 
function® for the operator Z and is important in the 
solution of differential equations with forcing functions. 

Comparison of«the coefficients of lower derivatives 
results in the further identities 


d 
a;(t)K(T, t) ie Hats t) ati Paice t), 


pS rd eon Bite (30) 


d 
ao(t) K(T, t) a ee ay 


The foregoing equations are the basis for a recurrence 
relation that permits the members of the fundamental 
set to be expressed in terms of the single function 
K(T, t), 


9) 


d 
qs en OK(T, DI. 


n—1 
Ua Dales al (31) 
w=) 
Hence, knowledge of the function K (¢, t’) is adequate to 
obtain the general solution to (19). 
The Green’s function K(f, ¢’) and the impulse re- 
sponse or weighting function R(t, ¢’) are equal through- 
out the interval of definition of the latter; that is, 


EP hE ed 


k(t, ’) = 0 ee (32) 


The description of linear control systems necessarily 
contains more complex equations than linear homo- 
geneous differential equations. One complicating factor 
is that most central systems utilize variations in a forc- 
ing function to achieve control. Hence, the preceding 
equations, which are a basis for the prediction of the 
transient component of the solution to a differential 
equation, must be augmented by formulas that include 
the effects of the steady-state component of the solu- 
tion, a component that arises from the presence of a 
forcing term. 


9K. S. Miller, “Properties of impulse responses and Green’s 
functions,” IRE TRAns., vol. CT-2, pp. 26-31; March, 1955. 


November 


The desired formulas can be derived in the following 
way. The differential 


n 


d’xpR 
= Do) —F = wl 


1=0 


L(xr) (33) 
has the function K(T, #) as an integrating factor. Multi- 
plication of both sides of (33) by K(T, t) and integration 
between the limits ¢ and T yield, according to (28), 


E pcr, tert) = xe(T) — | KT, t’)a,@)dt', (34) 


t 


or 


n—1 T 
wr(T) = >. f(T, Der(t) + K(T, #');(t) dt’... (35) 

2=0 t 
Whenever the forcing function x;(¢) is known over the 
interval to T, (35) provides an exact prediction formula. 
Such information usually is not available, and only 
estimates of the value of x;(¢) over this interval can be 
obtained. In Section II of this paper, it is shown that 
for optimum statistical prediction a “best” estimate of 
xs(t) is zero. In virtually every case of statistical optimi- 
zation the filter response is obtained from a combined 
integrodifferential operation on the filter error signal; 
that is, the forcing function x;(t) is derived from a differ- 
ential operation on the variable xz(t). In such cases, the 
differential pie 


n 


a ae ae a 


i=0 


L(ar) = =a oo = M(xz) (36) 
where m<n can be used to describe the filter. If again 
the function K(T, t) is used as an integrating factor, the 
prediction 


en(P) = DfT, )x2() 


+] K(T, | > b(t 


t i=0 


— | ae (37) 


can be derived. The integral in (31) can be simplified, 
if the interchange of summation and integration is al- 
lowable by successive integration by parts. This process 
is tedious and not repeated here. The results of the in- 
tegration by parts are expressible in the form 


aa(T) = Dir, tap (t) — Ss gi(T, Axe (t) 
“RMT, t')aen(t!) dt! (38) 
where 
EE = ye ys oe ()K(T, #’)] 
gi(T, t) = = Cuf(T, 1) 1=0,++-,m—1. (39) 


j=n—m+1 


£957 


The definition of unilateral Green’s function usually is 
extended to include K*(T, t’) as the impulse response of 
(36). Prediction often involves neglecting the integral 
term in (37), a process equivalent to the assumption that 
x(t’) is zero for t’ between ¢ and T. The remaining 
terms in (37) can be recombined by the aid of (39) into 
the form 


m—1 


we(T) = Dar, Nee) — ST Cof{T Da 
i=0 j=n—m+7 
a 
i Dslr. ier (t) — »> y CiOfAT, tee (0) 
= F jar, net) (40) 
where 
x(t) = xed), 7=0,---,n—m-—1 
= m2) - bo Cure (0 
t=n—m,---,n—1, (Af) 


Frequently, the quantities x,(#) can be identified with 
the outputs of the energy storage elements in the filter. 
Unfortunately, the quantities C;;(¢), which are com- 
binations of the functions a;(¢) and 6,(t), are difficult to 
determine, except for rudimentary systems. However, 
the weighting functions f;(T, ¢) depend only on the 
operator L, which is related to the signal shaping in an 
optimum filter and often can be found by standard 
methods. 

If K*(t, t’) is identified with the forward path of the 
optimum filter, then y(é, t’) is equal to K*(¢, t’) and 
(38) can be written 
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i=0 


— 3 g(r, t)xr™. 


t=0 


(42) 


This expression is a mathematical identity for all values 
of r=t. Thus, in particular, if 7 =t+-a, 


st | m—1 
Xrp(t) — Dy ft + a, t)xr™ — ne gilt + a, tar”, (43) 


7=0 1=0 
This expression provides a means for extrapolating over 
any nonnegative prediction interval a. A convenient 
form for a is a triangular wave of high frequency. Use of 
such a prediction interval would permit oscilloscope dis- 
play of continuous estimates of the future values of a 
- function. 
The synthesis procedure for the optimum extrapolator 
~ now can be summarized in the following steps: 
1) From the separable impulse response of the for- 
ward path of the optimum filter y(¢, 4), determine 
the differential equation. 
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2) Find the impulse response K(r, ¢) of the operator 
On XR. 

3) Find the functions f;(7, #) and g;,(r, ¢) in terms of 
K(r, t) and the coefficients of the differential oper- 
ator. 

4) Form (42) for the extrapolated response. 
Relationships derived in this section are not restricted 
to application in optimum systems alone but may be 
employed whenever analytic extrapolation of any func- 
tion generated by an integrodifferential operation is re- 
quired. One such application is found in final value con- 
trol systems. ?!° 


IV. 


The significance of the foregoing development may 
be elucidated by a thorough discussion of the physical 
meaning of the blocks in the diagram of Fig. 1. For the 
purposes of theoretical investigations, substantial sig- 
nificance can be attached to the mathematical results, 
which comprise a method for the solution to the predic- 
tion integral (7) when only the solution to the filter in- 
tegral (6) is known. However, a more complete under- 
standing of the operation of the filter and the extrap- 
olator enhances the usefulness of the method in prac- 
tical applications. For such a discussion, it is convenient 
to assume that the quantities «7g and x;y shown in Fig. 
2 are white noise and that the network shaping xx into 
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v(t,t)) 


SHAPING 
NETWORK 


Fig. 2—Input shaping when X, is white noise. 


the signal component of the input «gs has an impulse 
response v(t, #1). Thus, 


t 
xs(t) = f v(t, 1) xrg(t1)dty. (44) 
If 7 is a time greater than ¢, the present time, then 
xs(r) = if u(r, t1)Xrg(ti) dtr. (45) 


However, to evaluate this integral, values of «7s would 
be required over the interval between # and 7. Because 
these values are not known at time ¢, the best physically 
realizable approximation to (45) is 


xsp(T) 24) (7, ty) arg(t1)dty. (46) 


—~o 


10 C, W. Steeg and M. V. Mathews, “Final-value controller syn- 
thesis,” IRE Trans., vol. AC-2, pp. 6-16; February, 1957. 
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This portion of the integral (45) omitted in writing (46) 
is 


"5 


v(7, t1) Xrg(t1)dt, (47) 


vs(r) — sole) = f 


t 


which vanishes as ¢ approaches 7. The integral (47) also 
would vanish, making the approximation exact if xrs 
were identically zero over the interval between ¢ and r. 
Because xrg is assumed to be white noise with zero 
mean, the vanishing of the integral (47) is equivalent to 
replacing x;s by its mean value (zero) over the interval 
tSt,Sr. This conclusion is tantamount to the statement 
that the best prediction of white noise with zero mean 
is no prediction at all. 

If the shaping network were available for physical 
measurement of both xs and xyzs and their derivatives, 
the integral (46) could be replaced by its equivalent 
expression, analogous to (40), and the prediction made 
on the basis of these measurements. However, this pro- 
cedure must be abandoned in many physical problems 
such as an automatic radar-tracking system where the 
shaping network is a remotely located physical body, 
and only the output xg which might represent the posi- 
tion of the aircraft being tracked is available for meas- 
urement. Because such measurements inevitably are 
corrupted by noise, both filtering and prediction must 


be conducted on the basis of an error-minimization cri- ~ 


terion and associated optimization theory. 

The operation of the filter can be described in the 
following terms. The input to the filter is the noise-cor- 
rupted signal, and the response of the filter is the best 
approximation to the signal component of the input in 
the sense of least-mean-square error. If the noise com- 
ponent of the input is white, then the error signal x, is 
also white and has the same autocorrelation function 
as the noise. (This result is established in Appendix IT.) 
Consequently, the operation of the forward path of the 
filter, as shown in Fig. 3, is to recreate the signal com- 
ponent out of white noise. But this operation is pre- 
cisely the same as that of the shaping network which 
produces the signal originally. This reconstruction proc- 
ess can be regarded ascomposed of the three following 
steps: 1) a gain to change the amplitude of «xz into the 
amplitude of xrs, 2) a network to change xc into a sig- 
nal x4 with the optimum approximation to the time 
function of xzs, and 3) an exact model of the shaping 
network. Of these three steps, only the second involves 
any difficulty. Ideally, the second step would convert 
xc into xrg. The autocorrelation function of x¢ is the 
same as the autocorrelation function of xzs, but the 
corresponding time functions are different; therefore, 
this conversion cannot be performed exactly since 
neither time function is known. Hence, an error crite- 
rion, in this case, least-mean-square error, is used to 
minimize the difference between x;s and «4 as functions 
of time. The combination of the results of steps 2) and 
3) may be regarded as the optimized system model; that 
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Fig. 3—Alternate diagram for forward path of optimum filter. 


is, the model of the shaping network optimized with 
respect to minimization of the mean-square error. 

Because the forward path does represent a model of 
the shaping network, prediction can be made on the 
basis of signals physically existing in the filter, which is 
accessible for measurements, rather than on the actual 
system itself. This prediction is made by means of the 
relation 


Xrp(T) =i y(r, t1)“n(ti)dty (48) 


—o 


or its associated equivalent (41). Relation (48) is used 
instead of (46), which cannot be achieved, and the opera- 
tion of (41) is called extrapolation. 


V. CONCLUSION 


This paper presents a direct method for the time- 
domain synthesis of optimum extrapolators. The extrap- 
olators are optimum predictors that require merely the 
knowledge of design principles for optimum filters. To 
indicate the universal applicability of the synthesis pro- 
cedure, the method is used for the simplified design of 
predictors for two distinct situations: 1) prediction for a 
continuously varying instant that is always a fixed num- 
ber of seconds in the future, and 2) prediction for a 
specific fixed instant that is a continuously decreasing 
number of seconds in the future. The method is dis- 
cussed on the basis of uncorrelated signal and noise, 
both of which are derived from white-noise sources, al- 
though these restrictions can be removed easily by argu- 
ments analogous to the ones given by Bode and Shannon® 
in their discussion of prediction theory. Furthermore 
because the extrapolator design is independent of the 
error criterion, the synthesis procedure can be used for 
the simplified design of other systems such as predictors 
that use finite-memory filters.!! 

The hypothesis might be advanced that the time-de- 
pendent functions used for extrapolation are deter- 
mined solely from the signal characteristics. This con- 
jecture is true for the cases proved here and appears to 
hold under less restrictive conditions. If the extrap- 
olator functions are found by the method described, 
the effects of approximations inevitably introduced in 
the construction of optimum systems can be assessed 


nL. A. Zadeh and J. R. Ragazzini, “An extension of Wiener’s 
spent of prediction,” J. Appl. Phys., vol. 21, pp. 645-655; July, 
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readily. Although a completely general solution to the 
problem of optimum prediction still is lacking, the tech- 
niques for the design of optimum extrapolators may 
provide considerable insight into many previously un- 
solved cases and may aid in establishing rules for further 
understanding of optimum processes. The discussion is 
clarified by the following example. 


Example 1 
The configuration for a typical prediction system is 
obtained by assigning the following values: 
yis—ts(ti, ti + 7) = ors—rs(7) = 7Ad(z) 
yvn(i, 4 + 7) = onw(r) = 2ad(r7) 


o(t, #1) = (t — te), (49) 


This impulse response corresponds to a signal-shaping 
network described by 


tg + 2B%s3 + 6x3 = x78. 


The autocorrelation function for the signal component 
of the input is 


(50) 


1 
yss(t2, t) = 1A a. eAlt4l(g|[¢ — t2| +1). (51) 


Through the use of (66), the impulse response of the 
forward path of the optimum filter is obtained readily as 


1 
sg el eFC) [3(¢ — 12) +1] = way(t, ts), > te. (52) 


By means of (66), the solution to (7) has been circum- 
vented through choice of the unity-feedback configura- 
tion. The differential equation corresponding to this im- 
pulse response is 


PAR eee 
ée + 2Bxe + Pxr = — (=) (éz + 28xz), — (53) 
a \ 4p? 


a relation that can be derived by various techniques 


from (52). It is significant to note that the differential 
operator on the left-hand member of (53) is identical 
with the differential operator of (50). The integral (63) 
for g(t),!2 the closed-loop impulse response, can be 
written 


ft] 
eh) = 9h) — fo = Weleda. 8) 
0 
The solution to this equation, for 
cv fant 
yh) = ae 6t1(8t, 4-(1), (55) 


~ can be found directly, but the easiest method for solving 


this Volterra integral equation of the second kind is to 


a 


_ write the corresponding differential 


a Petes the filter has constant coefficients, the perrulec response 
is written as a function of a single variable. - 
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Afi 
ae SS (tr + 2Baxr). (56) 
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A 
As 
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A re-examination of (51) indicates that the forward 
path consists of the three parts indicated in the text: 1) 
a gain A/a, 2) a matching network (1/48?) (s+26), and 
3) a model of the signal shaping network. 

The functions required for the prediction operation, 
as indicated by (38), are 


fo(r, t) = e®-9[1 + B(r — 1] 
ita) = er SO G7) 


go(r, t) = Car ea (Te— “be (57) 


4aB? 
Verification of these functions is obtained by partial 
differentiation of (38) with respect to ¢ while 7 is held 
constant. From (43) the extrapolator for prediction at 
any nonnegative number: of seconds q@ in the future is 


A 
xrp = € E +d += Ba)te a ve (58) 


The extrapolator for prediction at a specific fixed future 
instant Tis 


tep = € B(P-o \(r — ioe + [1+ 6(T — dee 


rie dee Dah E60) 


The use of (58) with (53) yields the optimum pre- 
dictor transfer function 


; €*8[(1 + a8)s + B(2 + a8) | 
» (60) 


XRP i, Aa 
csannch ( : )s + B+ 
4a6? 2a? 


which agrees with results obtained by conventional 
methods. 

An analog-computer setup for this problem is shown 
in Fig. 4, where 8 has been allowed to approach zero in 
order to simplify equipment requirements. Although 
the filter shown is not strictly an optimum, the optimum 
is broad and the chosen configuration permits the merit 
of the method to be assessed more quickly. 

In Fig. 5, the typical results of an analog-computer 
study are shown for prediction at the fixed instant 
T =20 seconds. The function x7s(¢) consists of an impulse 
in the presence of white noise. The signal component of 
the input xs is shown on the first oscillograph record. 
The problem is to determine the value of «g at the time 
T =20 seconds when the filter input x7(t) is corrupted by 
white noise. The output of the optimum filter x(t) is 
shown in the second record where a slight displacement 
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Fig. 5—Analog computer results. 
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or delay in the response can be detected with respect to 
the signal component of the input. The bottom record 
shows the output of the optimum extrapolator. Al- 
though the peak-to-peak amplitude of the noise is ap- 
proximately 2.5 times the amplitude of the signal, 
xrp(t) rapidly assumes the correct mean value and the 
effects of the noise are nearly negligible over the last 3 
seconds. 


APPENDIX | 


In Section II, it was stated that the integral (16) 
must be neglected in determining future input values. 
This conclusion is given greater validity in the follow- 
ing argument. 

The definition (10) of the optimum filter allows (13) 
to be written as the difference of two integrals: 


x(t) = fv to) xr (te) dty 


-{ f y(t, ti) g(t, t2)ar(te)dtedty. (61) 


An interchange in the order of integration on the second 
integral permits the integrals to be recombined as one 
relation: 


val) = [y(t w)arltdat 


oe“ 


-{ | f y(t, ta) g(t, in) v(t 
=|] Ec te) -f y(t, tgttr, tn)dt arta, 


2 


(62) 
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Because this expression is the same as (10), the identity 
between impulse responses can be determined 


t 


g(t, te) = y(t, te) -f y(t, ti) g(t, te)dtr, 


t 


(63) 


an equation with a well-known counterpart in the theory 
of time-invariant systems.® When both sides of (61) are 
multiplied by x;(¢3) and averaged, 


t 
yrr(ts, ) = f y(t, to) yrr(ts, t2) dle 


—o 


t t2 
-f f y(t, to) g(te, t1)yrr(ts, t1)dtidte 


t t 
an oll, t)rnlts ta)dla = f y(t, te)yss(ts, te)dte 


=e —o 


= i y(t, tz)yww(ts, t2)dle. (64) 


The transition between the first and second equalities 
follows from the use of (7). Because the filter is optimum 


yrr(ts, t) = yss(ts, 2), boots (65) 


and 


t 
yss(ts, t) = it y(t, te)ynn(ts, t2)dle, t > ts (60) 


This equation frequently leads to an easy interpretation 
of the optimum system, particularly when the noise 
component of the input is white. This integral equation 
frequently is much easier to solve than (7), and the con- 
figuration of Fig. 1 was chosen in anticipation of the re- 
sult. 

If (14) is hypothesized to describe the predictor re- 
sponse, then, by use of (65), 


t 
y1—RP (to, t) = if y(t =f Qa, ti) yre(te, ty) dty 


—co 


t 
*. f y(t + a, tr)yrr(te, t1)dty 


00) 


t 
es f y(t + a, t1)yrr(te, t1)dty 


—3<O 


t 
-F i y(t + a, ti) ynw(te, t1) dtr. (67) 


When ¢ is replaced by t+a in (66) and comparison is 
made with (67), 


y1—rp(t2, t) = yss(te, t+ a) 
tra 
= f 96+ a tywalte,Addt (68) 
t 


When the noise component of the input is white, the 
latter integral vanishes, 


y1—rp(te, t) = yss(to, t+ a). (69) 


Loh 


Hence, by extending (63) the impulse response 
hit, te) = y(t + Qa, to) 


t 
— [vite melt tat, > by (70) 
tg 


can be introduced into (11) which becomes, by reason of 
(12), 


t 


Xxrp(t) ae i y(t + Qa, ti) x(t) dt, = f y(t -h Qa, t1) Xr (t1)dty 


— —oo 


t 
= || y(t + a, t)xn(ts)dty. 


a, 


(71) 


Eq. (71) provides a completely general relation be- 
tween the response of the optimum predictor and the 
error signal of the optimum filter in terms of the impulse 
response of the forward path of the optimum filter, 
whenever the noise component of the input is white. 


APPENDIX II 


In Section IV the statement is made that, if the noise 
component of the input is white, the error signal in the 
optimum filter has the same autocorrelation function. 
The proof of this contention is fundamental to the log- 
ical discussion of the validity of the extrapolationa 
technique. 

The autocorrelation function of the error signal is 


yer(ti, te) = xe(ti)xr(te) 
[xr(ts) — xe(tz) |[xr(t2) — xe(te) |. 


If g(t, ts) is the impulse response of the optimum filter 
defined by 


(72) 


t 
all) = fal, te)aa(ls dls (73) 
then g(t, ¢3) satisfies the integral equation for t>h 
t 
yss(ty) = [gtr 4d, (74) 


As a result of the use of (73) in (72), then for >h, 
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ty 
Veet, to) = yrr(ta, t1) =|} g(t1, ts) Xr(ts)xr(to) dts 


—o 


- fee ts) "r(t1)%1(ts) dls 


—eO 


te f g(t, ts) f ar(t3) Xr (ts) g(to, ty) dt3dt, 


ll 


te 
yir(ti, te) — if g(te, ts)yrr(ti, ts) dts 


—o 


= f Hes ts) Eze ts) 


t2 
rs { g(to, ta)yir(ts, tate dt. (75) 


08 


Because the signal and noise are assumed to be uncor- 
related, 


yir(ti, 4) = yss(t, t) + yuw(tr, 4), (76) 
and (74) can be written, for t>h, 
t 

yall ) = rally) = [g(t fret tate, (77) 


—o 


Thus, the right-hand side of (75) becomes, for t2>h, 


ty 
ya (li to) = ywvlly ) — [gt trwwlly todle. (78) 


If the autocorrelation function of the noise is 
ywn (ti, 41 + 7) = byn(7) = a8(7), (79) 


then, for r>0. 


1 
yer(h, 4 +7) = dnn(7) — f g(ti, ts)onn(li + 7 — ts)dts 


=O 


= onn(r) — ag(h, t: + 7). (80) 


However, the last term is zero because the impulse re- 
sponse g(t, t1) is zero for 4, >¢. Hence, 


venti, t+ 7) = den(t) = byn(z). (81) 
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